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Abstract:  This  study  characterizes  water  levels  in  the  Davis  Pond  floating 
marsh  created  by  the  diversion  of  fresh  water  from  the  Mississippi  River. 
The  model  was  validated  to  observed  field  data  collected  from  November 
2003  to  January  2004.  After  model  validation,  eight  alternatives  were 
tested  to  determine  the  benefits  of  extending  the  diversion  canal  (Alter¬ 
natives  1-3),  increasing  the  size  of  the  cuts  through  the  gabion  weir 
(Alternatives  4-6),  and  creating  breaches  in  the  Cypress  Lumber  Canal 
(Alternatives  7-8).  These  eight  initial  alternatives  were  analyzed  and  used 
to  create  four  additional  alternatives  (Alternatives  9-12)  consisting  of  the 
most  beneficial  aspects  of  each  of  the  initial  alternatives.  The  final  four 
alternatives  were  tested  to  determine  their  expected  benefits  to  the  system. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  Introduction 

The  Davis  Pond  freshwater  diversion  project  is  a  salinity-control  structure 
located  in  St.  Charles  Parish,  on  the  west  bank  of  the  Mississippi  River, 
two  miles  below  Luling,  LA.  Natural  and  man-made  changes  had  reduced 
freshwater,  sediment,  and  nutrient  inputs  from  the  Mississippi  River  to 
the  surrounding  estuarine-marsh  areas.  As  a  result,  saltwater  intruded 
further  inland  into  the  fresh,  intermediate,  and  brackish  water  zones, 
resulting  in  loss  of  habitat,  erosion,  and  deteriorating  water  quality.  The 
primary  objective  of  the  Davis  Pond  project  was  to  combat  the  saltwater 
intrusion  by  diverting  freshwater  from  the  Mississippi  River  into  the 
adjacent  estuarine  areas. 

The  Davis  Pond  project  was  dedicated  as  the  world’s  largest  freshwater 
diversion  in  March  2002.  Its  estimated  total  cost  was  $106  million 
(New  Orleans  District  2004).  The  project  was  designed  to  divert  up  to 
10,650  cubic  feet  per  second  (cfs)  of  freshwater  from  the  Mississippi  River, 
channeling  it  through  a  structure  via  an  outflow  canal,  into  a  9,200  acre 
ponding  area,  then  out  into  Lake  Cataouatche,  and  finally  into  the  remain¬ 
der  of  the  Barataria  Basin  (New  Orleans  District  2004).  The  ponding  area 
is  enclosed  on  the  north,  east,  and  west  by  guide  levees  with  design  eleva¬ 
tions  ranging  between  6.63  and  3.63  ft  NAVD881,  and  by  a  gabion  rock 
weir  on  the  southern  edge,  along  the  shoreline  of  Lake  Cataouatche.  At 
optimum  production,  the  Davis  Pond  diversion  would  preserve  33,000 
acres  of  wetlands  and  benefit  an  additional  777,000  acres  of  marshes  and 
bays  (Figure  1)  (New  Orleans  District  2004).  The  cost  benefits  is  estimated 
at  $15  million  for  fish  and  wildlife  and  an  additional  $300,000  for  recrea¬ 
tion  per  year  (New  Orleans  District  2004).  This  project  will  make  the 
Barataria  estuary  a  more  prolific  producer  of  oysters,  shrimp,  crab,  and 
fish,  as  well  as  a  major  habitat  for  migratory  waterfowl,  fur-bearing  ani¬ 
mals  and  alligators  (New  Orleans  District  2004).  The  increased  harvest 
(seafood  and  wildlife)  created  by  this  project  will  generate  significant 
economic  benefits  for  the  state  of  Louisiana. 


1  NAVD88:  North  American  Vertical  Datum  of  1988. 
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Figure  1.  Barataria  Basin  area  benefitting  from  the  Davis  Pond  diversion 
(New  Orleans  District  2004). 


When  normal  operations  of  the  Davis  Pond  freshwater  diversion  began, 
several  problems  were  identified  in  the  ponding  area.  The  first  and  most 
significant  was  that  stages  in  the  ponding  area  rose  much  higher  than 
anticipated.  This  issue  occasionally  resulted  in  undesired  flooding  when 
levee  overtopping  occurred.  The  problem  was  caused  in  part  by  a  previ¬ 
ously  unidentified  ridge  behind  the  gabion  rock  weir.  Also,  that  weir  was 
designed  with  a  specific  elevation  drop  to  compensate  for  expected  settling 
after  construction.  Actual  settling  of  the  weir  turned  out  to  be  significantly 
less  than  expected,  however,  resulting  in  less  flow  over  the  weir.  Another 
factor  limiting  flow  was  that  the  east  and  west  guide  levees  subsided  more 
than  expected,  and  were  overtopped  by  the  higher  stages  in  the  ponding 
area.  The  U.S.  Army  Engineer  District,  New  Orleans  (MVN)  cut  several 
notches  into  the  gabion  weir  in  an  attempt  to  lower  the  water  levels  and 
increase  the  amount  of  water  that  can  be  diverted.  While  that  action 
helped,  the  stages  in  the  ponding  area  were  still  too  high.  MVN  has  pro¬ 
posed  a  design  to  fix  these  problems  that  includes  widening  and  deepening 
the  existing  channels  through  the  rock  weir,  cutting  additional  channels, 
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extending  the  diversion  canal  farther  south,  cutting  breaches  into  the 
Cypress  Lumber  Canal,  and  raising  the  guide  levees  to  design  grades. 

MVN  requested  that  the  U.S.  Army  Engineer  Research  and  Development 
Center  (ERDC)  perform  a  two-dimensional  numerical  modeling  study  of 
the  ponding  and  surrounding  areas  to  assess  the  effectiveness  of  the 
proposed  designs. 
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2  Model  Development 

Description  of  site 

The  study  area,  shown  in  Figure  2,  is  located  in  southern  Louisiana,  south¬ 
west  of  New  Orleans.  Davis  Pond  is  located  primarily  in  St.  Charles  Parish 
with  a  small  area  residing  in  Jefferson  Parish.  The  Davis  Pond  freshwater 
diversion  is  intended  to  divert  flows  of  up  to  10,650  cfs  from  the  Missis¬ 
sippi  River  into  the  9,200  acre  Davis  Pond  marsh  area,  which  in  turn 
drains  into  Lake  Cataouatche.  The  Davis  Pond  area  is  covered  with  a  sig¬ 
nificant  amount  of  floating  marsh,  making  accurate  modeling  of  the 
hydrodynamics  of  the  area  very  challenging. 


Figure  2.  Outline  of  the  model  domain. 
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Field  data  monitoring 

Two  U.S.  Geological  Survey  (USGS)  gages  (Highway  90  and  Lake 
Cataouatche,  shown  in  Figure  9)  were  used  for  the  inflow  and  tidal  bound¬ 
ary  conditions.  MVN  also  installed  numerous  gages,  shown  in  Figure  3, 
throughout  the  Davis  Pond  marsh  area.  These  gages  were  vital  to  the  vali¬ 
dation  of  the  water  surface  elevations  in  the  model.  The  MVN  gages  were 
in  operation  from  November  30,  2003  to  January  13,  2004. 

It  was  determined  that  the  vertical  datums  for  the  observed  water  levels 
were  not  in  agreement.  Upon  further  investigation,  it  was  determined  that 
there  were  errors  in  the  vertical  datums.  Analysis  of  the  datums  is  pre¬ 
sented  in  Appendix  B,  along  with  the  method  used  to  determine  the  most 
likely  vertical  datums.  The  raw  observed  data  are  plotted  in  Figure  4,  with 
the  adjusted  observed  data  plotted  in  Figure  5.  These  figures  are  intended 
to  show  the  spread  in  water  surface  elevations  due  to  vertical  datum  issues 
for  the  raw  data  and  the  adjusted  data,  especially  around  hour  800. 


Figure  3.  Locations  of  water  surface  elevation  gages  used  to  validate  the  model. 
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Figure  4.  Raw  observed  data. 


Adjusted  Data 
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Figure  5.  Adjusted  observed  data. 
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Model  mesh  development 

The  numerical  model  was  developed  using  data  derived  from  several 
sources.  The  areas  north  of  the  gabion  rock  weir  were  defined  using 
LIDAR  data  and  topographic  surveys  provided  by  MVN.  The  LIDAR  data 
covered  the  areas  that  are  normally  dry  land  (elevations  above  o  ft 
NAVD88).  The  topographic  surveys  covered  the  east  and  west  guide  levees 
and  the  gabion  weir.  Bathymetric  surveys  covered  the  nine  channel  cuts  in 
the  gabion  weir,  which  extended  approximately  500  ft  into  the  marsh, 

Lake  Cataouatche  near  the  channel  cuts,  and  several  survey  cross-sections 
that  traversed  the  project  area.  Insufficient  surveys  were  available  to  do  a 
proper  interpolation  across  the  upper  and  lower  ponding  areas.  From  the 
available  cross-sectional  transects,  it  was  observed  that  the  ponding  area 
had  an  average  elevation  of  approximately  -1  ft  NAVD88.  Therefore,  this 
approximate  elevation  was  applied  to  the  entire  upper  and  lower  ponding 
areas.  The  horizontal  coordinate  system  used  in  the  model  was  State  Plane 
83,  South  Louisiana  zone  (ft),  with  the  vertical  datum  of  NAVD88  (ft).  The 
overall  model  domain  and  mesh  contours  are  presented  in  Figure  6. 


Figure  6.  Contour  plot  of  the  model  domain. 
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This  study  was  performed  using  the  TABS-MD  numerical  modeling  system 
(see  Appendix  C).  RMA2,  in  particular,  was  used.  The  specified  model 
coefficients  used  in  the  study  were  based  on  the  material  types  for  each 
element.  The  required  coefficients  specified  in  RMA2  include  an  eddy  vis¬ 
cosity  turbulent  mixing  coefficient,  a  bottom  friction  (Manning’s  n)  and 
the  coefficients  for  the  marsh  porosity  characteristics.  RMA2  version  5.0 
consists  of  a  new  layering  option  that  allows  for  the  specification  of  all 
required  coefficients  independently.  (Each  coefficient  has  its  own  layer.). 

Marsh  porosity  specification 

The  model  uses  a  technique  for  dealing  with  wetland  topography  and 
intertidal  marsh,  via  a  porosity,  to  allow  for  a  more  accurate  and  realistic 
representation  of  the  hydrodynamics  in  such  complicated  systems.  The 
marsh  porosity  was  specified  as  a  function  of  water  depth  such  that  no 
elevation  would  result  in  a  negative  depth  during  the  simulation.  A  more 
in-depth  explanation  of  marsh  porosity  and  the  modeling  of  floatant 
marshes  using  marsh  porosity  is  provided  in  Appendix  A.  The  specified 
marsh  porosity  coefficients  are  based  on  material  type.  Figure  7  shows  the 
marsh  porosity  materials  for  the  Davis  Pond  mesh.  Table  1  has  the  actual 
values  used  in  the  model  runs. 


MILES, 
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Figure  7.  Marsh  porosity  material  types. 
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Table  1.  Marsh  porosity  values. 


Material  Type 

AC1 

AC2 

AC3 

AC4 

1 

3.0 

2 

0.02 

0 

2 

3.0 

2 

0.02 

0 

3 

3.0 

2 

0.02 

0 

4 

3.0 

2 

0.02 

0 

5 

3.0 

2 

0.02 

0 

6 

3.0 

2 

0.02 

0 

7 

3.0 

2 

0.02 

0 

8 

3.0 

2 

0.02 

0 

9 

3.0 

2 

0.02 

0 

10 

4.0 

2 

0.02 

0 

11 

5.0 

2 

0.02 

0 

12 

6.0 

2 

0.02 

0 

13 

7.0 

2 

0.02 

0 

14 

8.0 

2 

0.02 

0 

15 

9.0 

2 

0.02 

0 

16 

10.0 

2 

0.02 

0 

17 

11.0 

2 

0.02 

0 

Note:  Marsh  porosity  values  used  in  the  numerical  model  (See  Appendix  A  for  an  explanation  of  the 
meaning  of  AC1,  AC2,  AC3,  and  AC4). 


Roughness  specification 

This  application  of  the  model  code  also  utilized  bottom  roughness  as  a 
function  of  the  water  depth.  The  mathematical  form  of  the  dependence  of 
the  Manning’s  friction  coefficient  with  depth  is 


n=—+nve  d/d° 
da  v 


(1) 


where: 

n0  =  scaling  friction  factor  for  depth  dependence 
d  =  water  depth 

a  =  exponent  on  depth  dependence 
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nv  =  scaling  factor  for  exponential  decay  dependence  (vegetative 
effects) 

do  =  reference  depth  for  exponential  decay. 

This  form  was  used  in  the  model  with  two  different  sets  of  parameters:  one 
in  the  channel  areas  and  a  separate  set  in  the  overbank  areas.  Figure  8 
shows  the  material  types  used  to  specify  the  frictional  values  in  the  model, 
with  Table  2  showing  the  actual  values  used  in  the  model  runs. 


Materials  Legend 

B  material  01 
material  02 
material  03 
material  04 
material  05 
material  06 
|  material  07 
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Figure  8.  Frictional  material  types. 


Table  2.  Frictional  values. 


Material  Type 

n0 

do 

nv 

a 

1 

0.06 

1.00 

0.10 

0.05 

2 

0.06 

1.25 

0.40 

0.20 

3 

0.06 

1.00 

0.40 

0.20 

4 

0.12 

2.00 

0.82 

0.41 

5 

0.12 

2.50 

0.82 

0.41 

6 

0.06 

1.25 

0.40 

0.20 

7 

0.12 

2.25 

0.82 

0.41 

ERDC/CHL  TR-08-11 


11 


Eddy  viscosity  specification 

The  eddy  viscosity  was  specified  using  the  Peclet  method.  The  equation  for 
eddy  viscosity  using  the  Peclet  method  is 


p  P  udx 
E 


(2) 


where: 

P  =  Peclet  number 
p  =  fluid  density 
u  =  average  elemental  velocity 
dx  =  length  of  element  in  streamwise  direction 
E  =  eddy  viscosity. 

A  global  value  of  15  was  used  for  the  Peclet  number  everywhere  in  the 
model  domain. 

Boundary  conditions  development 

The  boundary  conditions  for  the  model  consisted  of  a  tidal  boundary  at  the 
lower  end  of  Lake  Cataouatche  and  an  inflow  at  the  Davis  Pond  freshwater 
diversion.  Both  of  these  boundaries  were  updated  at  the  start  of  every  half- 
hour  time  step.  The  data  used  at  these  boundaries  were  obtained  from 
USGS  gage  measurements.  The  inflow  applied  was  the  observed  flow  at  a 
USGS  gage  located  in  the  diversion  canal  near  Highway  90.  The  Lake  Cata¬ 
ouatche  tidal  gage  was  located  on  the  eastern  shore  of  Lake  Cataouatche. 
The  locations  of  these  gages  are  shown  in  Figure  9.  A  plot  of  the  inflow  and 
tidal  boundary  conditions  applied  to  the  model  is  shown  in  Figure  10.  The 
Lake  Cataouatche  gage  data  plotted  in  Figure  10  were  obtained  by  shifting 
the  raw  Lake  Cataouatche  data  up  by  0.9  ft.  (See  Appendix  B  for  an  in- 
depth  explanation  of  the  reasoning  behind  this  shift.)  A  Hamming  filter 
was  also  applied  to  the  data  to  remove  any  noise  in  the  observed  signal. 
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Figure  9.  Location  of  USGS  gages. 


A  comparison  was  made  to  ensure  the  validity  of  applying  the  Lake  Cata- 
ouatche  observed  gage  values  at  the  southern  boundary  of  the  model 
domain.  The  plot  in  Figure  n  is  a  comparison  plot  of  the  Lake  Cataouatche 
observed  data  (applied  at  the  boundary)  and  the  model  results  at  the  loca¬ 
tion  of  the  Lake  Cataouatche  USGS  gage.  The  plot  shows  little  tidal  varia¬ 
tion  between  the  boundary  and  the  USGS  gage  location.  Therefore,  it  was 
assumed  that  no  significant  error  was  introduced  by  applying  the  Lake 
Cataouatche  gage  data  at  the  model  boundary.  The  boundary  conditions 
for  the  alternative  runs  consisted  of  a  steady-state  inflow  with  a  constant 
water  level  in  Lake  Cataouatche  of  1.3  ft  NAVD88  and  the  design  inflow  of 
10,650  cfs. 
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Boundary  Conditions 


—  Lake  Cataouatche  — Highway90  Discharge 
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Figure  10.  Plot  of  the  boundary  conditions. 


Figure  11.  Comparison  to  ensure  no  significant  effects  due  to  applying  the  Lake  Cataouatche 

gage  data  at  the  tidal  boundary. 
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Hydrodynamic  model  validation 

Comparisons  between  the  observed  gage  measurements  and  the  model 
results  are  shown  in  Figures  12-22.  The  pink  lines  are  the  adjusted  field 
data.  The  grey  haze  above  and  below  the  field  data  line  represents  the 
uncertainty  (±0.2  ft)  in  the  vertical  datums  (see  Appendix  B).  The  purple 
lines  are  the  model  results.  The  model  results  are  in  good  agreement  with 
the  observed  gage  data,  with  the  model  results  rarely  extending  beyond  the 
uncertainty  bounds.  The  gage  locations  were  shown  previously  in  Figure  3. 

The  comparisons  for  the  Highway  90  gage  and  Gage  23  (Figures  12  and  16, 
respectively)  are  the  only  ones  that  are  not  consistently  inside  the  error 
bounds.  It  is  believed  that  the  Highway  90  gage’s  close  proximity  to  the 
model  inflow  boundary,  along  with  the  crudeness  of  the  outflow  canal’s 
geometry  (a  trapezoidal  channel  approximated  as  a  square  channel),  is  the 
reason  for  this  unfavorable  comparison.  Because  Gage  20  (Figure  13) 
compares  very  well,  it  is  assumed  that  this  error  at  Highway  90  does  not 
propagate  very  far  into  the  system  and,  therefore,  is  of  little  concern  to  the 
overall  behavior  of  the  system.  Gage  23  is  located  close  to  the  gabion  weir, 
where  the  water  level  has  a  greater  slope.  Therefore,  any  slight  uncertainty 
in  the  gage  location  can  have  a  significant  effect  in  the  measured  water 
levels.  Also,  because  the  gage  is  near  to  the  gabion  weir,  it  is  possible  that 
supercritical  flow  is  occurring.  RMA2  is  unable  to  model  supercritical  flow, 
so  it  would  create  a  higher  water  level  to  compensate.  Since  Gages  27  and 
28  (Figures  17  and  18,  respectively)  compare  well,  it  is  believed  that  if 
there  is  an  error  in  the  model  results,  it  is  localized  at  the  gabion  weir  and 
does  not  affect  the  results  farther  into  the  system. 
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Validation  Comparison 


- Gage  20  Observed  - Gage  20  Model  Results 


Figure  13.  Plot  of  the  validation  comparison  at  Gage  20. 
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Figure  15.  Plot  of  the  validation  comparison  at  Gage  22. 
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Figure  17.  Plot  of  the  validation  comparison  at  Gage  27. 
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Validation  Comparison 


- Gage  29  Observed  - Gage  29  Model  Results 


Figure  19.  Plot  of  the  validation  comparison  at  Gage  29. 
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Figure  21.  Plot  of  the  validation  comparison  at  Gage  128. 
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Design  alternatives 

Alternatives  1-3  consisted  of  extending  the  diversion  channel  toward  the 
gabion  weir  using  different  channel  paths.  Alternatives  4-6  consisted  of 
modifications  to  the  outflow  channels  and  the  gabion  weir  to  allow  for 
more  flow  into  Lake  Cataouatche.  Alternatives  7-8  consisted  of  testing  the 
benefits  of  breaches  in  the  Cypress  Lumber  Canal.  The  first  eight  alterna¬ 
tives  were  used  to  test  the  effects  of  different  modifications  individually. 
From  the  results  of  the  eight  initial  alternatives,  four  additional  alterna¬ 
tives  were  created  to  determine  cumulative  effects.  Contour  plots  of  the 
bathymetry  of  the  different  alternatives  are  shown  in  Figures  23-34. 


ERDC/CHL  TR-08-11 


21 


Alternative  l  (Figure  23)  is  the  base/existing  case  (Figure  6)  with  the 
diversion  canal  extended  down  to  a  location  approximately  2,000  ft 
upstream  of  the  natural  ridge.  This  channel  configuration  was  chosen  in  a 
way  that  utilizes  existing  channels  as  much  as  possible.  The  channel  has  a 
bottom  width  of  120  ft  with  1  on  3  side  slopes.  The  initial  channel  invert  is 
-17  ft  NAVD88,  with  a  gradual  sloping  to  -18  ft  NAVD88  at  its  termination 
2,000  ft  upstream  of  the  natural  ridge. 
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Figure  23.  Contour  plots  of  Alternative  1. 
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Alternative  2  (Figure  24)  is  the  base  case  with  the  diversion  canal  extended 
down  to  a  location  approximately  4,000  ft  upstream  of  the  natural  ridge. 
This  channel  configuration  was  chosen  in  a  way  that  utilizes  existing 
channels  as  much  as  possible.  The  channel  has  a  bottom  width  of  120  ft 
with  1  on  3  side  slopes.  The  initial  channel  invert  is  -17  ft  NAVD88,  with  a 
gradual  sloping  to  -18  ft  NAVD88  at  its  termination  4,000  ft  upstream  of 
the  natural  ridge.  Alternative  2  is  Alternative  1  with  the  diversion  canal 
terminating  2,000  ft  farther  upstream  of  the  gabion  weir. 


Mesh  Elevation,  ft  NAVD88 

n3.00 
2.20 
1.40 
0.60 
-0.20 
|— |  -100 
-1.80 
-2.60 
H  -3.40 
—  4.20 
-5.00 


Mesh  Elevation,  ft  NAVD88 

n3.00 
2.20 
1.40 
0.60 
-0.20 
|— |  -100 
-1.80 
-2.60 
-3.40 

B-4.20 
-5.00 


Figure  24.  Contour  plots  of  Alternative  2. 


ERDC/CHL  TR-08-11 


23 


Alternative  3  (Figure  25)  is  the  base  case  with  the  diversion  canal  extended 
down  to  a  location  approximately  4,000  ft  upstream  of  the  natural  ridge. 
This  channel  configuration  is  a  straight  line  from  the  termination  of  the 
existing  diversion  channel  to  a  point  4,000  ft  upstream  of  the  natural 
ridge.  The  termination  point  for  Alternatives  2  and  3  are  identical,  with 
only  the  channel  paths  being  different. 
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Figure  25.  Contour  plots  of  Alternative  3. 
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Alternative  4  (Figure  26)  is  the  base  case  with  nine  85-ft-wide  cuts  in  the 
natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of  -1.5  ft 
NAVD88  (-6  ft  for  the  Cypress  Lumber  Canal  Cut).  The  approach  channels 
for  these  cuts  were  also  deepened  to  -1.5  ft  for  a  distance  of  approximately 
2,000  ft  upstream  of  the  gabion  weir.  The  side  slopes  for  the  channels 
were  1  on  3. 
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Figure  26.  Contour  plots  of  Alternative  4. 
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Alternative  5  (Figure  27)  is  the  base  case  with  nine  85-ft-wide  cuts  in  the 
natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  ft  for  the  Cypress  Lumber  Canal  cut).  The  approach 
channels  for  these  cuts  were  also  deepened  to  -3.0  ft  for  a  distance  of 
approximately  2,000  ft  upstream  of  the  gabion  weir.  The  side  slopes  for 
the  channels  were  1  on  3.  Alternative  5  is  Alternative  4  with  the  cuts  and 
channels  being  deepened  from  -1.5  to  -3.0  ft. 
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Figure  27.  Contour  plots  of  Alternative  5. 
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Alternative  6  (Figure  28)  is  the  base  case  with  nine  85-ft-wide  cuts  in  the 
natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  ft  for  the  Cypress  Lumber  Canal  cut).  Five  of  the 
approach  channels  for  these  cuts  were  deepened  to  -3.0  ft  for  a  distance  of 
approximately  2,000  ft  upstream  of  the  gabion  weir.  The  remaining  four 
channels  were  extended  even  farther  into  the  ponding  area.  The  side 
slopes  for  the  channels  were  1  on  3.  Alternative  6  is  Alternative  5  with 
three  of  the  approach  channels  extended  farther  into  the  ponding  area. 
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Figure  28.  Contour  plots  of  Alternative  6. 
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Alternative  7  (Figure  29)  is  the  base  case  with  one  breach  in  the  Cypress 
Lumber  Canal  Levee.  The  breach  was  85  ft  wide,  with  a  bottom  elevation 
of  -3  ft  NAVD88. 


MILES 
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Figure  29.  Contour  plots  of  Alternative  7. 
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Alternative  8  (Figure  30)  is  the  base  case  with  three  breaches  in  the 
Cypress  Lumber  Canal  Levee.  The  breaches  were  85  ft  wide,  with  a  bottom 
elevation  of  -3  ft  NAVD88.  Alternative  8  is  Alternative  7  with  two  addi¬ 
tional  breaches  in  the  Cypress  Lumber  Canal  Levee. 
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Figure  30.  Contour  plots  of  Alternative  8. 
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Alternative  9  (Figure  31)  is  the  base  case  with  nine  85-ft-wide  cuts  in  the 
natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  for  the  Cypress  Lumber  Canal  Cut)  with  1  on  3  side 
slopes.  The  approach  channels  for  these  cuts  were  also  deepened  to  -3.0  ft 
for  a  distance  of  approximately  2,000  ft  upstream  of  the  gabion  weir  for 
five  of  the  channels.  The  remaining  four  channels  were  extended  even 
farther  into  the  ponding  area.  There  are  also  three  breaches  in  the  Cypress 


Figure  31.  Contour  plots  of  Alternative  9. 


Lumber  Canal  Levee  in  the  same  locations  as  previously  identified  with 
Alternative  8.  This  alternative  also  includes  a  channel  connecting  the 
upper  and  lower  ponding  areas.  This  channel  has  a  -3.0  ft  NAVD88  eleva¬ 
tion  with  a  bottom  width  of  85  ft  and  1  on  3  side  slopes.  Alternative  9  is  a 
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combination  of  Alternative  6  (with  a  slightly  different  channel  alignment 
for  one  of  the  approach  channels)  and  Alternative  8  with  an  additional 
channel  cut  to  connect  the  upper  and  lower  ponding  areas. 

Alternative  to  (Figure  32)  is  the  base  case  with  eleven  85-ft-wide  cuts  in 
the  natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  for  the  Cypress  Lumber  Canal  Cut)  with  1  on  3  side 
slopes.  The  approach  channels  for  these  cuts  were  also  deepened  to  -3.0  ft 
for  a  distance  of  approximately  2,000  ft  upstream  of  the  gabion  weir  for 
six  of  the  channels.  The  remaining  five  channels  were  extended  even  far¬ 
ther  into  the  ponding  area.  There  are  three  breaches  in  the  Cypress 
Lumber  Canal  Levee  in  the  same  locations  as  previously  identified  with 


Figure  32.  Contour  plots  of  Alternative  10. 
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Alternative  8.  This  alternative  also  includes  a  channel  connecting  the 
upper  and  lower  ponding  areas.  That  channel  has  a  -3.0  ft  NAVD88 
elevation  with  a  bottom  width  of  85  ft  and  1  on  3  side  slopes.  Alternative 
10  is  a  combination  of  Alternative  6  (with  two  additional  cuts  and  slightly 
different  channel  alignments)  and  Alternative  8  with  an  additional 
channel  cut  to  connect  the  upper  and  lower  ponding  areas. 

Alternative  11  (Figure  33)  is  the  base  case  with  nine  85-ft-wide  cuts  in  the 
natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  for  the  Cypress  Lumber  Canal  Cut)  with  1  on  3  slopes. 
The  approach  channels  for  these  cuts  were  also  deepened  to  -3.0  ft  for 


Figure  33.  Contour  plots  of  Alternative  11. 
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approximately  2,000  ft  upstream  of  the  gabion  weir  for  five  of  the  chan¬ 
nels.  The  remaining  four  channels  were  extended  even  farther  into  the 
ponding  area.  There  are  three  breaches  in  the  Cypress  Lumber  Canal 
Levee  in  the  same  locations  as  previously  identified  with  Alternative  8. 
Alternative  11  also  includes  a  channel  connecting  the  upper  and  lower 
ponding  areas.  This  channel  has  a  -3.0  ft  NAVD88  elevation  with  a  bottom 
width  of  85  ft  and  1  on  3  side  slopes.  This  alternative  also  consists  of  two 
additional  channels  that  were  cut  from  the  locations  of  the  two  breaches  in 
the  Cypress  Lumber  Canal.  These  two  additional  channels  are  shown  in 
Figure  33.  They  have  bottom  elevations  of  -3  ft  NAVD88  with  side  slopes 
of  1  on  3.  Alternative  11  is  Alternative  9  with  the  two  additional  channels. 

Alternative  12  (Figure  34)  is  the  base  case  with  eleven  85-ft-wide  cuts  in 
the  natural  ridge  and  gabion  weir.  The  cuts  were  made  to  an  elevation  of 
-3.0  ft  NAVD88  (-6  for  the  Cypress  Lumber  Canal  Cut)  with  1  on  3  side 


Figure  34.  Contour  plots  of  Alternative  12. 
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slopes.  The  approach  channels  for  these  cuts  were  also  deepened  to  -3.0  ft 
for  a  distance  of  approximately  2,000  ft  upstream  of  the  gabion  weir  for 
six  of  the  channels.  The  remaining  five  channels  were  extended  even  far¬ 
ther  into  the  ponding  area.  There  are  three  breaches  in  the  Cypress 
Lumber  Canal  Levee  in  the  same  locations  as  previously  identified  with 
Alternative  8.  This  alternative  also  consists  of  two  additional  channels  that 
were  cut  from  the  locations  of  the  two  breaches  in  the  Cypress  Lumber 
Canal.  These  two  additional  channels  are  shown  in  Figure  34.  To  better 
capture  the  true  benefit  of  the  breaches  in  the  Cypress  Lumber  Canal,  the 
model  domain  was  extended  south  to  include  additional  channels  that 
branch  from  the  Cypress  Lumber  Canal.  These  channels  provide  addi¬ 
tional  connections  to  Lake  Cataouatche,  allowing  for  additional  outflow 
from  the  Davis  Pond  area.  The  channels  extending  south  from  the  Cypress 
Lumber  Canal  were  to  a  bottom  elevation  of  -3  ft  NAVD88  with  bottom 
widths  of  85  ft  and  side  slopes  of  1  on  3.  Alternative  12  is  similar  to 
Alternative  9,  but  without  the  channel  connecting  the  upper  and  lower 
ponding  area,  and  with  the  additional  domain  added  south  of  the  Cypress 
Lumber  Canal. 

Computational  environment 

The  hydrodynamic  model  runs  were  performed  on  the  ERDC  High  Per¬ 
formance  Computing  (HPC)  SGI  Origin  3000  (Ruby)  parallel-processing 
supercomputer.  The  base  conditions  model  contains  49,836  nodes  and 
17,893  elements.  The  model  was  executed  on  8  parallel  processors  and 
required  approximately  20  hours  of  computational  time  to  run  the  veri¬ 
fication  period  of  37  days  with  a  time  step  of  30  minutes.  The  production 
runs  required  a  shorter  duration  to  obtain  the  steady-state  solutions. 
Version  5.0  of  RMA2  was  used. 
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3  Hydrodynamic  Model  Results 

The  plot  in  Figure  35  shows  the  model  water  surface  profile  for  the  Davis 
Pond  study  area  for  Alternatives  1-8.  Figure  35  shows  Alternatives  4-6 
coming  closest  to  the  desired  profile,  while  Alternatives  1-3  and  7-8  pro¬ 
duce  insignificant  changes  to  the  water  levels  (Note:  some  of  the  profiles 
are  virtually  identical.)  The  desired  profile  was  provided  by  MVN.  There¬ 
fore,  it  is  reasonable  to  assume  that  the  primary  restriction  to  the  passage 
of  water  out  of  the  system  is  the  gabion  weir  along  the  southern  boundary 
of  the  Davis  Pond  marsh  area.  Additionally,  around  10,000  ft  upstream  of 
the  gabion  weir  there  is  a  sudden  increase  to  the  water  level  profile  (for 
Alternatives  4-6),  which  indicates  a  limit  in  the  channel  conveyance.  This 
result  corresponds  to  a  constriction  between  the  upper  and  lower  ponding 
areas  that  limits  the  benefits  obtained  from  Alternatives  4-6. 


These  initial  results  were  analyzed  to  create  Alternatives  9-12.  The  model 
results  for  these  three  additional  alternatives  are  shown  in  Figure  36  and 
Figure  37.  While  these  alternatives  result  in  significant  reductions  in  water 
surface  elevation,  the  water  levels  for  all  alternatives  are  still  significantly 
higher  than  the  desired  water  surface  elevation  profile.  Since 
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Alternative  12  was  deemed  the  most  desirable,  additional  runs  were  done 
with  inflows  of  7,500  and  5,000  cfs  for  this  case. 


Water  Surface  Profile  for  Davis  Pond  Marsh  Area 


- Desired  Profile 

- Base  Conditions 

Alternative  12  5000  cfs 

Alternative  12  7500  cfs 

Alternative  12  10650  cfs 

Figure  37.  Water  surface  profile  obtained  from  Alternative  12. 
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The  previous  plots  depict  profiles  along  a  specific  path  for  each  alternative. 
Figures  38-50  show  contour  plots  of  water  surface  elevations  for  the 
entire  ponding  area.  These  plots  provide  a  more  comprehensive  idea  of  the 
total  changes  experienced  by  the  system.  Difference  plots  of  the  base  con¬ 
ditions  and  Alternatives  9-12  are  shown  in  Figures  51-52.  These  figures 
show  the  benefits  that  would  be  obtained  by  the  implementation  of 
Alternatives  9-12. 


Figure  38.  Water  surface  elevation  for  the  base  conditions  in  feet  NAVD88. 
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Figure  39.  Water  surface  elevation  for  Alternative  1  in  feet  NAVD88. 


Figure  40.  Water  surface  elevation  for  Alternative  2  in  feet  NAVD88. 
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Figure  41.  Water  surface  elevation  for  Alternative  3  in  feet  NAVD88. 


Figure  42.  Water  surface  elevation  for  Alternative  4  in  feet  NAVD88. 
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Figure  43.  Water  surface  elevation  for  Alternative  5  in  feet  NAVD88. 


Figure  44.  Water  surface  elevation  for  Alternative  6  in  feet  NAVD88. 
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Steady  State  Solution  with  Q=1 0,650  cfs 


Figure  45.  Water  surface  elevation  for  Alternative  7  in  feet  NAVD88. 


Figure  46.  Water  surface  elevation  for  Alternative  8  in  feet  NAVD88. 
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Figure  47.  Water  surface  elevation  for  Alternative  9  in  feet  NAVD88. 


Steady  State  Solution  with  Q=1 0,650  cfs 


Figure  48.  Water  surface  elevation  for  Alternative  10  in  feet  NAVD88. 
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Figure  49.  Water  surface  elevation  for  Alternative  11  in  feet  NAVD88. 


Steady  State  Solution  with  Q=1 0,650  cfs 


Figure  50.  Water  surface  elevation  for  Alternative  12  in  feet  NAVD88. 
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Base  Solution  -  Alternative  9  Solution 

Base  Peak  Water  Level  -  Alternative  1 0  Peak  Water  Level,  ft 


Base  Solution  -  Alternative  10  Solution 


Figure  51.  Contour  difference  plots  of  the  base  conditions  and  Alternatives  9-10. 
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Base  Peak  Water  Le\^l  -  Alternative  1 1  Peak  Water  Level,  ft 
■  2.0 
5  1-8 
—  1.6 
—  1.4 
—1 1.2 
1.0 
0.8 
H  0.6 
0.4 
0.2 
0.0 


Base  Solution  -  Alternative  11  Solution 


Base  Solution  -  Alternative  12  Solution 


Figure  52.  Contour  difference  plots  of  the  base  conditions  and  Alternatives  11-12. 
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4  Conclusions 

From  the  initial  set  of  alternatives,  it  was  determined  using  the  hydro- 
dynamic  model  that  the  primary  restriction  to  the  passage  of  water 
through  the  system,  and  the  concomitant  excessive  water  levels,  was  the 
gabion  weir  in  conjunction  with  the  ridge  just  north  of  the  weir  located 
along  the  southern  boundary  of  the  Davis  Pond  marsh.  The  ridge  and  rock 
weir  restricted  the  flow  of  water  from  the  marsh  area  into  Lake 
Cataouatche,  creating  perched  water  levels  in  the  Davis  Pond  marsh. 
Alternatives  4-6  increased  the  size  of  the  cuts  in  the  ridge  and  weir,  allow¬ 
ing  more  water  to  pass  into  Lake  Cataouatche.  It  was  noted  that  these 
increased  connections  to  Lake  Cataouatche  revealed  a  channel  conveyance 
restriction  between  the  upper  and  lower  ponding  areas  (see  the  water 
surface  elevation  profile  plot,  Figure  35,  for  Alternatives  4-6  at  10,000  ft 
upstream  of  Lake  Cataouatche).  Alternatives  7  and  8  produced  little 
improvement,  but  that  result  could  be  due  the  model  domain  limitations. 
Alternative  12  was  tested  with  areas  added  to  the  model  domain  to  include 
channels  to  the  south  in  order  to  allow  for  additional  means  of  conveying 
water  from  the  Cypress  Lumber  Canal  into  Lake  Cataouatche.  These 
additional  areas  allowing  water  to  pass  out  of  the  ponding  area  made 
Alternative  12  the  most  beneficial  alternative  in  terms  of  water  level. 

No  alternative  produced  results  at  or  below  the  desired  water  surface  ele¬ 
vation  profile.  The  most  extreme  alternatives  produced  water  surface  ele¬ 
vation  profiles  that  were  a  foot  or  higher  than  the  desired  water  surface 
elevation  profiles.  Alternative  12,  with  an  inflow  of  5,000  cfs,  produced 
results  that  were  still  over  0.25  ft  above  the  desired  profile.  Additional 
lowering  of  the  discharge  could  produce  a  water  surface  profile  at  or  below 
the  desired  profile,  but  this  lower  inflow  may  not  meet  the  habitat  restor¬ 
ation  goals  for  the  entire  Barataria  Basin. 

From  these  results,  it  can  be  concluded  that  the  desired  water  surface 
elevation  profile  is  an  unrealistic  goal  for  anything  short  of  drastic 
modifications  to  the  system  or  drastic  lowering  of  the  desired  inflow. 
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Appendix  A:  Description  of  Marsh  Porosity 
and  Its  Application  to  Floatant  Marsh 

General  overview 

The  marsh  porosity  method  was  developed  (Roig  1995)  as  a  transition 
method  between  wet  and  dry  states  within  RMA2.  The  partial  motivation 
was  to  eliminate  the  severe  shocks  that  occur  with  conventional  wetting 
and  drying,  which  must  rather  arbitrarily  define  elements  as  either  wet  or 
dry.  Combine  the  binary  wet/dry  state  decision  with  numerical  oscilla¬ 
tions,  and  the  model  often  becomes  unstable  while  attempting  to  resolve 
wet/dry  boundaries  within  the  nonlinear  iterations  for  a  time  step  for  the 
conventional  wetting  and  drying  logic. 

The  basic  philosophy  of  the  marsh  porosity  method  is  that  there  is  a  grad¬ 
ual  and  continuous  variation  between  the  wet  and  dry  states  that  is  both 
more  realistic  and  more  stable  numerically.  The  implementation  of  the 
method  involves  the  definition  of  a  fractional  wetted  area  as  a  function  of 
water  surface  elevation.  The  implementation  of  the  method  within  RMA2 
is  limited  to  a  four-phase  variation  (see  Figure  Ai)  for  the  case  of  a  falling 
water  surface  elevation.  In  the  first  phase,  at  an  extremely  high  relative 
water  level,  the  node  is  viewed  as  completely  wet,  with  a  fractional  wetted 
area  of  1.0  (above  an  elevation  of  4  ft  in  Figure  Ai).  As  the  water  level  falls, 
the  node  enters  the  second  transition  phase,  where  the  fractional  wetted 
area  decreases  linearly  with  falling  water  level  (between  elevations  1.2  and 
3.0  in  Figure  Ai).  The  third  phase  is  the  portion  of  the  vertical  distribution 
where  the  fractional  wetted  area  takes  on  a  minimum  value  (AC3  in 
Figure  Ai  for  elevations  between  0.0  and  1.2).  The  final  phase  would  be 
the  dry  phase  when  the  water  level  falls  below  0.0.  In  the  actual  imple¬ 
mentation  in  RMA2,  the  drying  logic  will  dry  the  node  when  the  local 
water  depth  falls  below  the  drying  threshold  (DSETD  on  the  DE  card). 

The  marsh  porosity  method  is  implemented  within  RMA2  on  a  nodal 
basis.  This  has  been  the  source  of  considerable  confusion  to  users  develop¬ 
ing  an  understanding  of  the  technique.  The  fractional  wetted  area  is  more 
easily  conceptualized  on  an  elemental  basis  but,  in  fact,  it  is  actually  repre¬ 
sentative  of  the  local  area  associated  with  each  node.  Mathematically, 
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Marsh  Porosity 


Figure  Al.  Basic  marsh  porosity  specification  in  RMA2. 


the  influence  of  the  marsh  porosity  of  each  node  is  distributed  over  each 
element  based  on  the  finite  element  nodal  basis  functions  of  the  element. 
This  means  that  the  closer  one  is  to  the  node,  the  more  effective  its  nodal 
marsh  porosity  value  is,  and  the  farther  away,  the  less  influence  it  has. 

Ultimately  the  decision  of  whether  an  element  is  wet  or  dry  comes  down  to 
the  status  of  the  wet  or  dry  state  of  each  of  the  nodes  in  the  element.  In  the 
conventional  wetting  and  drying  approach,  without  marsh  porosity,  an 
element  is  made  dry  if  a  single  node  within  the  element  becomes  dry. 
When  marsh  porosity  is  invoked,  the  element  becomes  dry  only  after  every 
node  in  the  element  is  “dry”.  This  is  logical  in  that  if  it  were  the  same  as 
the  conventional  method  of  wetting  and  drying,  then  it  would  not  provide 
the  desired  transition.  By  deferring  the  drying  process  until  all  nodes  are 
dry,  the  shock  is  removed.  In  the  process  of  deferring  drying,  the  remain¬ 
ing  wet  nodes  within  the  element  are  passing  less  and  less  water  due  to  the 
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restricting  effect  of  reducing  the  fractional  wetted  area.  The  flow-restrict¬ 
ing  effect  is  accentuated  when  the  bottom  roughness  is  prescribed  as  a 
function  of  water  depth. 

Mathematical  implementation 

The  mathematics  of  the  implementation  are  rather  elegant.  The  basic 
governing  equations  for  RMA2  are  summarized  below,  as  extracted  from 
RMA2  users’  guide  (USACE  2003). 


The  generalized  computer  program  RMA2  solves  the  depth-integrated 
equations  of  fluid  mass  and  momentum  conservation  in  two  horizontal 
directions.  The  forms  of  the  solved  equations  are 
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where: 

h  =  depth 

u,v  =  x  and  y  direction  velocities,  respectively 
x,y,t  =  Cartesian  coordinates  and  time 
p  =  density  of  fluid 
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8  =  eddy  viscosity  coefficient, 

for  xx  =  normal  direction  on  x-axis  surface; 
yy  =  normal  direction  on  y-axis  surface; 
xy  and  yx  =  shear  direction  on  each  surface 
g  =  acceleration  due  to  gravity 
a  =  elevation  of  bottom 
n  =  Manning’s  n  value 

1.486  =  conversion  from  SI  (metric)  to  non-SI  units 
£  =  empirical  wind  shear  coefficient 
Va  =  wind  speed 
=  wind  direction 


co  =  rate  of  earth’s  angular  rotation 
O  =  local  latitude. 


The  application  of  the  marsh  porosity  method  is  accomplished  by  defining 
an  effective  marsh  porosity  depth.  This  is  defined  as  the  integral  of  the 
fractional  wetted  area  from  the  relative  zero  elevation  in  Figure  Ai  (bottom 
of  the  offset  channelized  area)  to  the  water  surface  elevation: 

hn  =  f  Kdz  (A2) 

Jo 


where  K  is  the  fractional  wetted  area  presented  in  Figure  Ai.  The  specifica¬ 
tion  of  K  is  as  follows  for  the  special  case  that  h  >  ztop: 


=  AC  3 


=  AC3  + 


for  z<zbot 

for  Ztoi<z<z„ 
for  ztop  <z  <h 


(A3) 


For  this  special  case  (h  >  ztop )  the  effective  depth  can  be  evaluated. 
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In  the  application  of  this  to  the  governing  equations,  there  are  two  steps 
that  should  be  mentioned.  First,  the  depth  h  defined  above  is  a  modified 
depth  from  the  depth  displayed  within  the  Surface-water  Modeling  System 
(SMS)  (Brigham  Young  University  [BYU]  2002)  for  the  nodal  elevation. 
The  first  thing  that  is  done  when  marsh  porosity  is  used  is  to  “offset”  the 
bottom  elevation  by  the  distance  ACi  (h  =  Hsms  -  ACi).  The  next  step  is  to 
develop  an  effective  depth  factor  based  on  the  integrated  value  above. 


If  the  governing  equations  above  are  multiplied  through  by  o,  the  equa¬ 
tions  are  converted  to  the  same  equations  as  above,  but  with  ha  substituted 
for  h.  The  wind  stress  and  friction  terms  will  be  appropriately  reduced  by 
a,  which  varies  between  o  and  1,  corresponding  to  the  reduced  surface  area 
over  which  these  terms  act.  When  the  nonlinear  equations  are  solved  via 
Newton-Raphson  iteration,  the  Jacobian  matrix  has  additional  terms 
dh 

involving  — - .  When  the  water  surface  is  well  above  ztop  these  corrections 
dh 

are  not  significant.  Equation  A4  shows  that  the  marsh  porosity  results  in  a 
reduction  in  the  actual  water  depth.  If  the  coefficients  ACi,  AC2  and  AC3 
are  zero,  then  zbot  =  ztop  =  A0,  the  elevation  from  the  SMS  geometry  file, 
and  ha  =  h. 

The  usefulness  of  the  marsh  porosity  method  has  been  extending  beyond 
its  original  intent.  It  was  recognized  that  the  vertical  variation  in  the 
marsh  porosity  factor  is  directly  analogous  to  the  statistical  distribution  of 
bottom  elevation.  With  this  analogy  complex,  sub-mesh  scale  geometric 
variations  can  be  represented  on  various  spatial  scales  using  marsh  por¬ 
osity.  This  approach  was  first  applied  to  wetlands  on  the  coast  of  Louisiana 
on  spatial  scales  of  miles  to  represent  canals  and  levees  along  with  natural 
wetland  topography  (Letter  1993).  It  has  been  routinely  applied  to  other 
estuarine  systems. 

The  use  of  marsh  porosity  as  a  statistical  representation  of  sub-grid  scale 
features  makes  the  use  of  conventional  wetting  and  drying  less  of  a  neces¬ 
sity.  Because  the  method  now  represents  the  statistical  variation,  physi¬ 
cally  removing  the  element  from  the  computational  domain  has  less  mean¬ 
ing.  This  leads  to  an  additional  alternative  in  the  application  of  marsh 
porosity,  where  the  offset  elevation  (ACi)  is  set  such  that  the  statistical 
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minimum  elevation  is  always  wet  and  no  wetting  and  drying  is  performed. 
This  approach  is  very  stable  numerically. 

Floatant  marsh 

The  study  area,  Davis  Pond,  is  covered  extensively  by  floatant  marsh, 
where  water  is  able  to  flow  through  and  underneath  a  floating  mat  of  vege¬ 
tation.  The  effect  of  this  on  the  routing  of  flow  through  the  system  was 
evaluated  by  considering  the  analogy  to  marsh  porosity. 

The  floatant  marsh  is  considered  in  the  context  of  marsh  porosity  in  two 
ways.  First  the  analogy  of  the  presence  of  the  floatant  marsh  on  a  uniform 
horizontal  mineral  bottom  is  considered.  Then  the  presence  of  the  floatant 
marsh  in  addition  to  the  variable  bottom  elevation  with  a  description  via 
mash  porosity  is  evaluated. 

Assume  for  consideration  a  thickness  of  the  floating  marsh,  D,  of  which  a 
thickness  d  is  submerged.  The  water  depth  is  h  above  the  mineral  bottom. 
The  porosity  of  the  floating  mat  is  Pm  as  shown  in  Figure  A2. 


Figure  A2.  Schematic  of  floatant  marsh  parameters. 


The  effective  flow  area  per  unit  width  will  be 

h~d  +  Pmd 


Expressing  this  area  as  a  fractional  area  of  the  full  depth  area  gives  a 
“marsh  porosity  factor”  of 
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K  =  ±-d{±-Pm)  /  h  forh>d 
=  Pm  forh<d 


(A5) 


This  factor  is  plotted  in  Figure  A3  for  d  =  1.0  and  Pm  =  0.15.  The  curve  is 
nonlinear  for  values  of  h  >  d.  In  an  attempt  to  linearize  the  distribution  for 
application  in  RMA2  we  can  assume  that  after  the  depth  gets  larger  than 
the  submerged  portion  of  the  floating  mat  (h>  d),  the  flow  is  restricted  to 
the  area  below  the  mat  (h  -  d).  The  mat  may  have  some  small  movement, 
but  that  flow  will  be  offset  by  the  reduction  in  the  velocity  profile  adjacent 
to  the  bottom  of  the  mat  (see  Figure  A2).  We  will  therefore  assume  that 
there  is  no  contribution  of  flow  from  the  water  within  the  mat.  This  will 
simplify  the  marsh  porosity  distribution  to 

K  =  1  —  d  /  h  for  h  >  d 

(A6) 

=Pm  for  h  <  d 


Figure  A3.  The  marsh  porosity  factor  for  example  d  =  1.0,  Pm  =  0.15. 
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For  this  simplification  to  work,  we  will  need  to  recognize  the  further 
restriction  that  the  value  of  K  must  be  greater  than  Pm.  This  is  also  non¬ 
linear,  and  is  presented  in  Figure  A2.  Note  that  for  this  simplification,  the 
marsh  porosity  reaches  the  minimum  value  of  Pm  at  a  greater  water  depth 
than  for  Equation  Ai.  This  effect  is  more  severe  for  the  case  of  Pm  =  0.4 
with  still  d  =  1.0,  as  presented  in  Figure  A4.  This  approach  will  therefore 
not  be  reasonable. 


Figure  A4.  Marsh  porosity  representation  for  of  =  1.0,  Pm  =  0.4. 


So  to  assume  a  linear  form  of  the  marsh  porosity  factor,  we  assert  the  fol¬ 
lowing  form: 


K  =  Pm  +  a (h  -  d)(  1  -  Pm )  still  for  h  >  d  (A7) 

In  order  to  completely  define  this  distribution  we  must  select  which  is  in 
fact  the  inverse  of  AC2  in  the  RMA2  input.  For  the  case  presented  in 
Figure  A3,  we  set  a  =  1/3. 
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The  selection  of  an  appropriate  value  of  a  could  come  from  a  minimization 
of  the  deviation  from  the  theoretical  nonlinear  curve.  It  could  be  estimated 
in  one  of  several  ways.  First,  the  minimization  integral  can  be  defined  as 

J^(k„-kb)dz  (A8) 

where  the  kn0n  is  the  nonlinear  relation  in  Equation  At  and  kun  is  the  linear 
relationship  in  Equation  A2.  If  we  assume  we  know  the  maximum  water 
depth  c/max,  we  can  set  the  integral  to  zero.  The  difficulty  in  evaluating  this 
integral  is  that  the  limits  (up  to  dmax)  will  likely  involve  the  extension  of 
Equation  A2  to  depths  for  which  k  >  1.  To  accurately  evaluate  the  integral 
requires  that  the  integral  be  split  into  two  pieces: 

/!'“  (*™,  -  hn  )dz  =  P”  (km  -  k,in  )dz  +  /*“  (km  - 1  )dz  (A9) 

Jd  Jd  Jztop 


where  ztop  =  the  top  elevation  of  the  ultimate  linear  transition 
(Ao  +  AC2/2).  Unfortunately,  ztop  is  determined  by  a,  which  we  are 
attempting  to  solve  for.  Therefore,  we  would  need  to  minimize  the  inte¬ 
grals  in  an  iterative  process. 


— (1  -P 

h[ 


-(Pm+a(h-d)(l-Pm))\dh 


— (l-p  ' 

h[  m/ 


-l\dh  =  0 


(A10) 


(h-d) 


dh  =  0 


These  integrals  can  be  evaluated  to  solve  for  a  as 


a  = 


2a -pj 

(■ ztoP-d ) 


1 


(All) 
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Once  a  is  estimated,  then  ztop  can  be  corrected  as  ztop  =  d  +  l/a,  as  shown 
in  Table  Ai.  The  estimate  for  ztop  is  then  best  taken  as  the  average  of  the 
previous  value  and  the  corrected  value.  For  our  example  above,  if  we  use 
d  =  l.o,  Pm  =  0.4  and  assumed  hmax  =  5,  using  an  initial  estimate  of  ztop  as  4 
we  iterate  to  a  value  of  0.409  for  a  and  ztop  as  3.44. 


Table  Al  Iterative  solution  of  Equation  All. 


Estimate  for  ztop 

Computed  value  of  a 

Corrected  value  of  ztop 

4 

0.556225 

2.797834 

3.398917 

0.394918 

3.532173 

3.465545 

0.416674 

3.399958 

3.432751 

0.406115 

3.462358 

3.447555 

0.410916 

3.433585 

3.44057 

0.408658 

3.447034 

3.442293 

0.409216 

3.443696 

3.442994 

0.409443 

3.44234 

3.442667 

0.409338 

3.442972 

There  is  an  alternative  way  to  optimize  the  linear  fit.  Expressing  the  depth 
as  a  function  of  k  for  the  transition  range  of  depth  gives 


d(l  ~PJ 


(i-9 


hlm  -  - 


(k-P,n) 

“(1  -PJ 


(A12) 


The  selection  of  an  appropriate  value  of  a  could  come  from  a  minimization 
of  the  deviation  from  the  theoretical  nonlinear  curve.  Now  the  minimiza¬ 
tion  integral  can  be  defined  as 

J^iKon-WJdk  (A13) 

•Mnin 


Notice  that  the  integral  is  limited  to  some  maximum  value  kmax.  This  is 
required  because  the  nonlinear  depth  estimate  has  a  singularity  at  k  =  1. 
Therefore,  if  we  assume  a  maximum  marsh  porosity  kmax  (corresponding  to 
the  value  at  dmax),  we  can  set  the  integral  to  zero  and  solve  for  a.  In  this 
case  we  do  not  have  to  break  the  integral  into  two  pieces. 
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X 


fcmax  |  Cl(l  ~  Pm  ) 


d  + 


(k-PJ 


a(l  -Pm)j 


(1  -*) 

This  integral  can  be  evaluated  and  solved  for  a  as 

(k  —  P  )2 

\  max  m ) 


dk  =  o 


(A  14) 


a  = 


2d(l-P„ 


1  — Pjln 


1 

-P 

m 

1- 

b 

max 

—  k  4-  P 

max  1  m 


(A15) 


The  best  way  to  approach  the  estimate  of  /cmax  is  to  use  equation  l  with  hmax 
for  h.  Using  the  example  from  above,  the  estimated  km ax  would  be  0.88. 
This  yields  a  value  of  a  as  0.395,  as  compared  with  0.409  from  the  method 
above. 

These  values  of  a  can  then  be  used  to  estimate  the  marsh  porosity  values 
for  use  in  RMA2  for  the  flat  mineral  bottom  case. 


AC1  =  d 


AC  2 
AC3 


1 

a(l-Pm) 

Pm 


(A  16) 


Table  A2  presents  calculations  with  varying  values  of  d,  Pm  and  hmaxfor  the 
flat  bottom  handling  of  the  floatant  marsh. 
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Table  A2  Floatant  marsh  on  flat  bottom.  Variation  of  AC2  computations 
for  variations  in  d,  Pm  and  h  max. 
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Combined  effect  of  floating  marsh  and  conventional  marsh  porosity 

The  normal  specification  of  marsh  porosity  will  complicate  the  incorpora¬ 
tion  of  the  effects  of  the  floating  marsh.  Therefore,  for  completeness  we 
will  evaluate  the  effects  of  the  floating  mat  in  combination  with  conven¬ 
tional  marsh  porosity  specification. 

There  are  several  cases  that  need  to  be  considered.  These  involve  the  rela¬ 
tive  floating  mat  thickness  compared  with  the  variation  parameter  (AC2 
on  the  DM  card)  for  conventional  marsh  porosity,  and  the  water  depth. 
The  goal  of  this  analysis  is  to  develop  the  functional  relationship  between 
the  effective  marsh  porosity  (i.e.,  including  floatant  marsh  effects)  and 
water  depth.  These  cases  can  be  grouped  by  where  the  water  surface  eleva¬ 
tion  lies  within  the  vertical  variation  of  the  marsh  porosity  (see  Figure  Ai). 
These  cases  are: 

1.  Case  1 :  the  water  surface  is  higher  than  ztop. 

2.  Case  2:  the  water  level;  falls  in  the  range  Zbot  <h<  ztop. 

3 .  Case  3 :  the  water  level  is  lower  than  Zbot. 

4.  Case  4 :  the  water  surface  is  below  zero  and  the  node  is  dry  (not 
considered). 

Each  case  has  sub-cases  based  on  how  deep  the  bottom  of  the  floating  mat 
reaches  within  the  marsh  porosity  variation.  These  cases  are  defined  in 
Table  A3  and  illustrated  in  Figure  A5. 


Table  A3.  Cases  for  floatant  marsh  in  conjunction  with  marsh  porosity  bottom. 


Case 

Water  Surface  Elevation 

Bottom  of  Floating  Mat 

la 

h  >  ztop 

h  -  d  >  ztop 

lb 

h  >  Ztop 

Zbot  <h  -d  <  ztop 

lc 

h  >  ztop 

0  <h  -d  <  Zbot 

Id 

h  >  Ztop 

h  -  d  =  0 

2a 

Zbot  <  Ztop 

Zbot  <h  -  d  <  ztop 

2b 

Zbot  <  Ztop 

0  <h  -d  <  Zbot 

2c 

Zbot  <11  <  Ztop 

h-d  =  0 

3a 

0  <h  <  Zbot 

0  <h  -  d  <  Zbot 

3b 

0  <h  <  Zbot 

h  -  d  =  0 
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Figure  A5.  Definition  of  sub-cases  for  combined  floating  marsh 
with  marsh  porosity  bottom  elevation  variation. 
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Case  la 

This  gives  the  marsh  porosity  factor  vertical  distribution  as 


K 


AC  3 


AC  3  4 


(1  —  AC3)(z- zb) 
AC  2 


1 

Pm 


for  z<zb 

for  zb<z<zt 

for  zt  <z  <{h-d) 
for  Qi-d)<z  <h 


The  integrated  composite  marsh  porosity  is  then 


w;^=r“+/: 

J>  h-d  n h 

„  dz+L/-dz 


AC  3T 


(l-AC3)(z-zJ 


AC  2 


dz 


=  AC3(zA  +  h-d-zt+Pmd+ 


(1-AC3)AC2 


For  this  case,  the  effective  value  of  K  is 


ldK 

h  dh 


1 


Case  lb 

When  the  water  level  drops  such  that  the  floating  marsh  mat  straddles  the 
upper  limit  of  the  elevation  zt,  but  is  higher  than  Zb  (see  Case  lb  in 
Figure  A5),  we  will  have 


K 


AC  3 


AC  3  + 


1  —  AC3)[z  —  zt 


AC  3 


AC  2 

1  —  AC3)[z  —  zb) 


AC  2 


for  z<zb 
for  zb<z<  h-d 

for  {h-d)<z<zt 


for  zt  <  z  <h 


The  composite  marsh  porosity  effective  depth  is  then 
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h°=fo  Kdz  =  Jo  AC3dz+l 

(l  —  AC3)(z  —  zb 


h-d 

zb 


AC  3  4 


(l-AC3)(z-zJ 


AC  2 


\dz 


+ 


AC  3  + 


r 

J  h-d 


f  b AC3dz+  f  ' 

Jo  J  zb 


AC  2 


AC  3  4 


Pmdz  +  £Pmdz 


(l-AC3)(z-zh) 


AC  2 


\dz 


J'zt 

h-d 


AC  3  + 


1  —  AC3)[z  —  zt 


AC  2 


1  -  P,n  )  dz  +  /‘/!Pmdz 


=  AC3  ■  zt  +  Pmd  —  AC2 


-^(1-AC3)  +  AC3 

2  v  ’ 


+hb  [  AC3  +  Pm(l-  AC3)]  -  (1  -  Pm )  (12  ^  K 


where  hb  =  h  -  d  -  Zb. 


For  this  case,  the  effective  composite  value  of  K  is 


K  - dh° 
Keff  ~  dh 


AC3+PJ1 


x  1  — Pm  1-AC3  , 

AC3)  -  - — - ’-hh 

AC  2  b 


This  equation  will  apply  as  long  as  the  floating  mat  straddles  the  upper 
marsh  porosity  elevation  zt,  regardless  of  the  thickness  of  the  mat.  If  the 
thickness  of  the  floating  marsh  mat  is  greater  than  AC2  we  will  obtain 
Case  lc. 


Case  lc 


=  AC3 

for 

zmin  <z<h-d 

=  AC3  Pm 

for 

h-d  <z  <  zbot 

=  P„,{aC3+(1^3)(z 

for 

Zbot  <~-  z  <~-  ztop 

=  Pm 

for 

ztoP  <z<h 

Integrating  to  define  the  effective  water  depth: 
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J'h  nh-d  rzhnt 

Kdz  =  /  AC 3  dz+  AC 3  Pm  dz  + 

0  Jo  Jh-d  m 

+ \aC3  +  (1~^3)(Z-Zi «)W  +f  P-dz 

U  zbot  -/1U  Z  U  ztop 


-  AC 3  (l -Pm)(h-d)  +  Pmztop  +  Pm 


1-AC3)AC2 


+  h  —  z. 


For  this  case,  the  effective  value  of  K  is 


Case  Id 


(711 

K«f=-^i=AC3+pJ1-ACi) 


=  AC 3  P 


for  0  <z  <zh 


=  Pm  |  AC3  +  \^\z  -  zbot)  1  for  zbot  <z<  ztop 


for  ztov  <z  <h 


Integrating  to  define  the  effective  water  depth: 


K  =  AC3Pmdz 

+ f*Pm  Uc3  +  +  f  Pmdz 

U  zbot  -/lUZ  U  zto p 


1-AC3)AC2 


+  h-(l-AC3)z. 


For  this  case,  the  effective  value  of  K  is 


K  _  dho  _  p 
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Case  2a 


K 


=  AC  3 


=  AC  3 


(1  —  AC3) 


z  — 


=  ^  MC3 


AC2 
(1-AC3) 
AC2 


f0r  Zmin  <  Z  <  Zbot 

for  zbot  <z  <h-d 
for  h  —  d<z<h 


Integrating  to  define  the  effective  water  depth: 


Xh  f*zhnt  r»h—d  \ 

j  Kdz  =JQ  AC 3  dz+J  j  AC 3  + 

+fi/fAC3+ 


h-d  I 

zbot 

(1  — AC3) 


(1-AC3) 


AC  2 


z-zbot)\dz 


AC  2 


)\dz 


=  AC3(h-d  +  Pmd ) 


(1-AC3)  r 


2  AC2 


{(!  -  pm  ){h  —  d-  zbot )  +Pm(h-  zbot  y 


For  this  case,  the  effective  value  of  K  is 


K  _  dh0 
Kff  ~  dh 


AC3  +  (\c23){/l  ~ ' d - Z“  +  P'"dJ 


Case  2b 


K= 


AC  3 
AC3  Pm 


|aC3  + 


(1-AC3) 

AC2 


M  zmin  <z  <h-d 
for  h-d  <  z  <  zhot 

for  zbot  <z  <h 


Integrating  to  define  the  effective  water  depth: 
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nh  nh-d  rzhnt 

h=  Kdz  =  /  AC3  dz  +  /  AC 3  Pdz 

°  Jo  Jo  Jh-d  m 

SLP~ 


+ 


AC3{h-(l-Pm)d}- 


(1-AC3) 

2AC2 


Pm{h-Zboty 


For  this  case,  the  effective  value  of  K  is 


K  _dK 

Kff  ~  dh 


AC3+P 


(1-AC3) 

AC2 


{' h-Zbot } 


Case  2c 


=  AC 3  Pm  for  zmin  <  z  <  zbot 

=  pm  |^C3  +  (1^4C23)(Z Zfeo4  f°r  Zbot  <Z<h 


Integrating  to  define  the  effective  water  depth: 


h°=S0Kdz=ffAC3P™dz 

+r  PiAC3+<kAm 

Jzbot  m  AC 2 


Z-Zbot)\dz 


=  pm\AC3h+{12^\h-zboty 


For  this  case,  the  effective  value  of  K  is 


v  _  dK 

Keff~~dh 


AC  3  4 


(1-AC3) 


AC  2 


{h-zbot} 


Case  3a 


K= 


AC  3 
AC3  Pm 


M  zmin  <  z  <h-d 
for  h  —  d<z<h 


Integrating  to  define  the  effective  water  depth: 
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ha  =  fh Kdz  =  fh  d  AC3  dz  +  AC3  Pmdz 

°  Jo  Jo  Jh-d  m 

=  AC3{h-(l-Pjd} 


For  this  case,  the  effective  value  of  K  is 


(>K 

dh 


AC  3 


Case  3b 


K=  AC  3  Pm  for  zmin  <z  <h 


Integrating  to  define  the  effective  water  depth: 


K  =  /;  Kdz  =  1C3  Pmdz  =  AC3  Pm  h 


For  this  case,  the  effective  value  of  K  is 


dh 

K*=^h=AC3P" 


The  equations  for  ha  and  Kejf  are  summarized  in  Table  A4  for  each  of  the 
cases.  At  the  bottom  of  the  table  are  the  equations  for  the  cases  without 
floatant  marsh.  These  can  be  confirmed  to  be  the  same  for  each  subcase  by 
setting  the  mat  porosity  to  a  value  of  1.  The  value  of  Ke/f  for  the  Case  1 
subcases  all  is  1  for  Pm  =  1. 


Table  A4.  Summary  of  combined  floatant  marsh  with  marsh  porosity  for  bottom. 


AC3(zt  )  +  h  —  d  —  zt  +Pmd+  (1  -  AC3)AC2  /  2 

Acm  +l^gr+  AC3W  -  ac^ac2Ia> 

2AC2  m  ‘  2AC2 


r,  w  ,  i  (l-AC3)AC2 

AC3[(1  - Pm)(h -d)  +  Pmztop ]  +  Pm  y- - -> - +  h-zt( 


1-AC3)AC2 

Pm  - J1 - +  ^  -  (1  ”  AC3)Z« 


AC3(h-d  +  Pmd)+ 


(1  —  AC3) 
2AC2 


{(!  -  pm  ){h  -  d  -  zbot  f  +  Pm  (h  -  zbot  y 


AC3[k  -  (1  -  P„  )d}  +  (12  y™  f„  (ft  - , 

F"\Ac3h+il^(h-ZbA 

AC3{h-(l-Pm)d} 

AC3  Pm  h 

n=1)  (l  —  AC3)AC2  /  2  +  /i  —  (1  —  AC3)ztop 

"=1)  ,  (1-AC3),, 

+  2  AC2  ( k  Zbot  ^ 

n=l)  AC3  h 


AC3  +  P  (1  —  AC3)— 


AC3+P  1-AC3) 


(1-PJ(1-AC3) 


Aca+^-^ih-d-^  +  p^} 

'4C3  +  ^C1AC23){'1~Z-'} 
P-[j4C3  +  (1AC23){'!~Z“} 


AC3  P, 


'4C3+(1AC23)|'1~Z“> 
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Summary 

The  previously  discussed  derivations  have  shown  that  floatant  marsh  can 
be  represented  using  marsh  porosity  if  extensive  information  on  the  char¬ 
acteristics  of  the  area  is  available  (i.e.,  marsh  thickness,  horizontal  cover¬ 
age,  and  frictional  resistance).  For  the  Davis  Pond  study,  none  of  these 
characteristics  was  known,  so  it  was  deemed  unnecessary  to  incorporate 
the  additional  capabilities  into  RMA2.  The  existing  configuration  is 
believed  to  produce  an  accurate  representation  of  the  behavior  of  the 
system.  Modification  of  the  code  to  achieve  full  representation  of  the 
floatant  marsh  would  require  input  quantities  for  the  unknown  marsh 
thickness,  coverage,  and  frictional  resistance.  Therefore,  the  RMA2  marsh 
porosity  was  not  modified  and  was  used  as-is  for  this  project. 
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Appendix  B:  Datum  Analysis 

Lake  Cataouatche  tidal  datum 

In  organizing  the  data  for  the  Davis  Pond  verification,  comparisons  were 
performed  on  the  Highway  90  and  Lake  Cataouatche  U.S.  Geological 
Survey  (USGS)  gages.  The  locations  of  these  gages  are  shown  in  Figure  Bi, 
with  the  observed  data  plotted  in  Figure  B2.  The  U.S.  Army  Engineer  Dis¬ 
trict,  New  Orleans  (MVN)  reported  to  the  Engineer  Research  and  Develop¬ 
ment  Center  (ERDC)  that  the  Highway  90  gage  was  returning  water  sur¬ 
face  elevations  that  were  high  by  0.132  ft.  This  shift  was  applied  to  the 
gage  data,  and  henceforth  any  discussion  of  the  Highway  90  water  surface 
elevation  data  incorporates  this  correction.  The  plot  in  Figure  B2  shows  a 
significant  head  difference  from  Highway  90  to  Lake  Cataouatche  even 
after  prolonged  periods  with  little  to  no  inflow  from  the  Davis  Pond  diver¬ 
sion,  a  situation  that  seems  unlikely.  It  was  decided  that  additional  investi¬ 
gation  should  be  performed  to  determine  if  this  reported  behavior  is  an 
accurate  representation  of  the  true  performance  of  the  system. 


Figure  Bl.  Location  of  the  Highway  90  and  Lake  Cataouatche  USGS  gages. 
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The  gage  analysis  was  performed  using  a  two-pronged  approach  with  com¬ 
parisons  being  made  with  other  USGS  gages  and  National  Oceanic  and 
Atmospheric  Administration  (NOAA)  gages.  Since  each  organization  has 
different  personnel  and  methods  for  determining  the  vertical  datum  of  its 
gages,  any  agreement  between  USGS  and  NOAA  gages  produces  a  high 
confidence  level  for  accuracy.  The  locations  of  these  USGS  and  NOAA 
gages  are  shown  in  Figure  B3.  The  data  for  the  USGS  gages  were  provided 
in  a  gage  datum  with  the  conversion  from  this  gage  datum  to  NAVD88  also 
being  provided  by  the  USGS.  The  data  for  the  NOAA  gages  were  obtained 
from  the  NOAA  website.  The  Port  Fourchon  and  Grand  Isle  historical  data 
were  downloaded  in  Mean  Sea  Level  (MSL).  A  conversion  from  MSL  to 
NAVD88  was  obtained  from  the  NOAA  website  for  the  Grand  Isle  gage. 
(Port  Fourchon  had  no  conversion  to  NAVD88.)  Due  to  their  close  prox¬ 
imity,  this  conversion  was  applied  to  the  tidal  data  for  both  gages. 
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Figure  B3.  Location  of  gages  used  to  investigate 
the  Lake  Cataouatche  vertical  datum  issues. 

The  Lake  Cataouatche  gage  data  are  plotted  with  the  other  USGS  gages, 
Lake  Salvador  and  Little  Lake,  in  Figure  B4;  and  with  the  NOAA  gages, 
Grand  Isle  and  Port  Fourchon,  in  Figure  B5.  Both  plots  are  in  agreement 
that  the  Lake  Cataouatche  measurements  are  low.  Short-term  variations  in 
tides  and  phase  differences  make  it  difficult  to  determine  an  appropriate 
shift  using  the  raw  tidal  data.  To  avoid  these  complications,  a  filter  was 
applied  to  the  gage  data  to  remove  the  high-frequency  signals  (Periods 
<  1,000  hr)  from  the  data.  This  left  only  the  extreme  low-frequency  sig¬ 
nals,  which  tend  to  be  a  representation  of  the  mean  water  level  of  the  tidal 
signal.  Plots  of  the  raw  data  and  the  low  frequency  signals  are  shown  in 
Figures  B6  (USGS  gages)  and  B7  (NOAA  gages).  It  should  be  noted  that 
there  were  gaps  in  the  Lake  Salvador  data.  The  method  used  to  fill  these 
gaps  to  provide  the  most  accurate  filtered  low-frequency  signal  is  dis¬ 
cussed  in  the  Lake  Salvador  Gage  Data  Analysis  section. 
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Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Raw  Data 


Figure  B4.  Plot  of  the  Lake  Cataouatche  gage  data 
with  the  data  from  other  USGS  gages. 


Figure  B5.  Plot  of  the  Lake  Cataouatche  gage  data 
with  the  data  from  the  NOM  gages. 
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Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Raw  Data  and 
Low  Frequency  Signal  (Periods  >  1,000  hrs) 
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Figure  B6.  Plot  of  the  raw  data  and  low-frequency  signals 
for  Lake  Cataouatche  and  USGS  gages. 


Comparison  of  Lake  Cataouatche  with  NOAA  Gages.  Raw  Data  and 
Low  Frequency  Signal  (Periods  >  1,000  hrs) 
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Figure  B7.  Plot  of  the  raw  data  and  low-frequency  signal 
for  Lake  Cataouatche  and  the  NOAA  gages. 
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From  these  low-frequency  signals,  an  upward  shift  of  0.9  ft  was  deter¬ 
mined.  By  applying  this  shift,  the  Lake  Cataouatche  data  lay  between  the 
two  NOAA  gage  and  between  the  two  additional  USGS  gages.  The  shifted 
low-frequency  signals  are  shown  in  Figures  B8  (USGS)  and  B9  (NOAA). 
While  this  shift  is  believed  to  be  the  most  likely  vertical  datum  level  for  the 
Lake  Cataouatche  gage,  it  is  by  no  means  exact.  Even  with  the  extensive 
analysis  performed  on  the  Lake  Cataouatche  data,  the  vertical  datum 
should  only  be  considered  accurate  to  within  ±0.2  ft.  This  0.2  ft  uncer¬ 
tainty  was  chosen  since  in  Figures  B8  and  B9,  all  filtered  data  signals 
seemed  to  be  within  ±0.2  ft  of  the  shifted  Lake  Cataouatche  filtered  signal. 
Due  to  this  uncertainty,  the  model  was  validated  using  the  0.9  ft  shift 
determined  from  this  analysis,  and  then  a  sensitivity  analysis  was  per¬ 
formed  to  determine  the  importance  of  the  0.2  ft  uncertainty  associated 
with  the  vertical  datum  for  the  Lake  Cataouatche  gage.  Plotted  in  Fig¬ 
ures  Bio  and  B11  are  the  filtered  and  unfiltered  signals  with  the  shifted 
Lake  Cataouatche  data  as  compared  with  the  other  USGS  and  NOAA 
gages.  It  should  be  noted  that  this  procedure  was  applied  to  more  current 
data  (2007)  and  it  showed  no  significant  difference  in  mean  water  level 
between  the  Lake  Cataouatche  gage  and  the  remaining  gages  used  in  this 
analysis.  Therefore,  it  is  assumed  that  the  vertical  datum  issue  addressed 
here  has  been  corrected. 

Plotted  in  Figure  B12  is  the  Highway  90  data  with  the  shifted  Lake 
Cataouatche  data.  There  is  no  longer  a  significant  head  difference  between 
the  Highway  90  and  Lake  Cataouatche  gages  during  long  periods  of  little 
to  no  inflow  from  the  Davis  Pond  diversion.  The  remaining  gages  in  the 
Davis  Pond  marsh  area  were  shifted  up  or  down  to  lay  at  the  same  vertical 
level  as  Highway  90  and  Lake  Cataouatche  during  long  periods  of  little  to 
no  inflow.  This  was  deemed  appropriate  since  there  should  be  no  high 
water  surface  elevations  in  between  two  low  water  surface  elevations  for 
this  system.  The  raw  and  adjusted  data  are  plotted  in  Figures  4  and  5, 
presented  previously  in  the  main  text. 

Lake  Salvador  gage  data  analysis 

The  raw  Lake  Salvador  gage  data  are  plotted  in  Figure  B13.  There  are  two 
extended  gaps  in  the  data  (hours  1609  to  2235  and  hours  2350  to  2821).  If 
there  were  no  knowledge  of  the  behavior  of  the  system,  the  usual  approxi¬ 
mation  would  be  to  enter  the  mean  of  the  data  for  these  gaps,  but  by  using 
the  Lake  Cataouatche  data,  a  better  approximation  can  be  achieved. 
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Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Low  Frequency 

Signal 

(Periods  >  1,000  hrs) 


Time,  hr  (0  is  Oct  1  0:00.  2003) 

Figure  B8.  Low-frequency  signal  for  USGS  gages  and  shifted  Lake  Cataouatche  data. 


Comparison  of  Lake  Cataouatche  with  NOAA  Gages,  Low  Frequency 

Signal 

(Periods  >  1,000  hrs) 


Figure  B9.  Low-frequency  signal  for  NOAA  gages  and  shifted 
and  unshifted  Lake  Cataouatche  data. 
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Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Raw  Data  and 
Low  Frequency  Signal  (Periods  >  1,000  hrs) 


Figure  BIO.  Raw  data  and  low-frequency  signal  for  USGS  gages 
and  shifted  Lake  Cataouatche  data. 


Comparison  of  Lake  Cataouatche  with  NOAA  Gages,  Raw  Data  and 
Low  Frequency  Signal  (Periods  >  1,000  hrs) 


Figure  Bll.  Raw  data  and  low-frequency  signal  of  NOAA  gages 
and  shifted  Lake  Cataouatche  data. 
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Comparison  of  Lake  Cataouatche  with  USGS  Highway  90  Gage 
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Figure  B12.  Plot  of  the  Highway  90  and  Lake  Cataouatche  gage  data. 


Figure  B13.  Lake  Salvador  raw  data. 
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Figure  B14  displays  a  plot  of  the  Lake  Salvador  data  and  the  Lake 
Cataouatche  data.  With  the  exception  of  the  vertical  datum  shift,  these  two 
signals  are  very  similar.  The  Lake  Cataouatche  signal  was  shifted  so  as  to 
best  match  the  Lake  Salvador  signal  (up  0.7  ft).  The  shifted  Lake 
Cataouatche  signal  was  used  from  hour  1609  to  hour  2821  in  place  of 
missing  Lake  Salvador  data.  Figure  B15  shows  the  raw  signal  and  the  new 
signal  using  the  Lake  Cataouatche  data.  To  perform  a  sensitivity  test  on 
the  effects  of  this  missing  data  to  our  low-frequency  signal,  a  filter  was 
applied  to  the  created  Lake  Salvador  data  set  discussed  above,  and  a  filter 
was  also  applied  to  the  Lake  Salvador  data  set  with  the  mean  value 
inserted  in  the  gaps.  These  two  filtered  data  sets  are  plotted  in  Figure  B16. 

While  there  is  a  definite  difference  in  the  low-frequency  signals,  the  differ¬ 
ence  appears  to  be  mainly  relegated  to  -300  hours  before  and  after  the 
interruption  in  data.  Both  filtered  signals  are  very  similar  for  hours  o  to 
1300  and  hours  3100  to  4500.  Therefore  the  values  inputted  for  the  data 
gaps  have  a  limited  effect  on  the  data  prior  to  hour  1300  and  after 
hour  3100. 


Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Raw  Data 


Figure  B14.  Lake  Salvador  and  Lake  Cataouatche  raw  data. 
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Comparison  of  Lake  Cataouatche  with  USGS  Gages,  Raw  Data 


Figure  B16.  Lake  Salvador  raw  signal  and  created  signal  with  filtered  signals. 
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Effects  of  frictional  uncertainty 

An  analysis  was  performed  to  determine  an  approximate  effect  of  the 
±0.2  ft  uncertainty  in  the  vertical  datum  used  to  validate  the  numerical 
model.  This  ±0.2  ft  uncertainty  could  have  affected  the  frictional  values 
obtained  from  the  model  validation.  Therefore,  the  model  frictional  values 
were  perturbed  to  create  water  surface  elevations  that  were  increased  and 
decreased  by  ±0.2  ft.  As  discussed  in  Chapter  1,  the  frictional  values  are 
separated  into  two  categories:  overbank  and  channel.  There  were  three 
options  for  obtaining  the  ±0.2  ft  perturbation  in  water  levels: 

1.  Increase/decrease  the  overbank  frictional  values. 

2 .  Increase /decrease  the  channel  friction  values . 

3 .  Increase/decrease  all  frictional  values . 

Options  1  and  2  are  assumed  to  produce  the  most  extreme  results.  Both 
possibilities  were  examined.  For  option  1,  the  base  conditions  were  run 
with  the  steady-state  plan  boundary  conditions  to  obtain  the  most  likely 
results.  The  overbank  frictional  values  were  perturbed  until  the  water  level 
was  raised  and  lowered  by  0.2  ft.  This  resulted  in  a  high  bound,  “most 
likely”  level,  and  a  lower  bound  for  the  base  conditions.  To  determine  the 
approximate  effect  on  the  alternative  water  surface  elevations,  Alternative 
10  was  run  with  each  of  the  three  sets  of  frictional  values.  Alternative  10 
was  chosen  because  it  was  one  of  the  more  extreme  alternatives.  This  same 
procedure  was  used  for  option  2  with  the  channel  frictional  values.  Plots  of 
the  water  surface  elevation  profiles  for  this  analysis  are  shown  in  Figures 
B17  and  B18. 

The  plots  in  Figures  B17  and  B18  show  the  importance  of  variation  in 
overbank  frictional  values  versus  channel  frictional  values.  It  took  little 
perturbation  with  the  overbank  frictional  values  to  create  the  0.2  ft  high 
and  low  bounds  (110  %  and  90  %  frictional  values,  respectively),  while 
both  channel  frictional  values  took  significantly  more  perturbation  (130% 
and  60%)  for  the  base  conditions.  From  this  analysis,  it  can  be  concluded 
that  the  overbank  frictional  values  play  a  larger  role  than  the  channel 
frictional  values  in  the  water  surface  elevations  for  the  base  conditions. 

But  because  the  channel  frictional  perturbation  has  the  larger  range  with 
Alternative  10,  the  channel  frictional  values  become  more  important  after 
the  modifications  have  been  made  to  the  system.  Therefore,  the  water  in 
the  system  has  gone  from  moving  through  the  overbank  areas  to  moving 
primarily  through  the  channels  in  the  system. 
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Modification  of  Overbank  Friction  to  Determine  the  Importance  of  the 
0.2  ft  Uncertainty  in  the  Vertical  Datums 


—  Base  Conditions  —Base  Conditions  with  Lower  Overbank  Friction 

—  Base  Conditions  with  Higher  Overbank  Friction  - Alternative  10 

- Alternative  10  with  Lower  Overbank  Friction  —Alternative  10  with  Higher  Overbank  Friction 


Figure  B17.  Plot  of  the  effect  of  the  0.2  ft  uncertainty  on  the  vertical  datums 
by  modifying  the  overbank  frictional  values. 


Modification  of  Channel  Friction  to  Determine  the  Importance  of  the 
0.2  ft  Uncertainty  in  the  Vertical  Datums 


- Base  Conditions  —Base  Conditions  with  Lower  Channel  Friction 

—  Base  Conditions  with  Higher  Channel  Friction  —Alternative  10 

—  Alternative  10  with  Lower  Channel  Friction  —  Alternative  10  with  Higher  Channel  Friction 


Figure  B18.  Plot  of  the  effects  of  the  0.2  ft  uncertainty  in  the  vertical  datums 
by  modifying  the  channel  frictional  values. 
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This  analysis  determined  that  the  planned  water  surface  elevation  profiles 
could  be  shifted  up  or  down  by  as  much  as  a  half-foot  in  Alternative  to 
based  on  the  uncertainty  of  the  vertical  datum  of  the  gage  data  used  to 
validate  the  numerical  model.  It  should  be  noted  that  the  0.5  ft  uncertainty 
for  the  alternative  water  surface  elevation  profiles  would  apply  to  the  most 
extreme  combination  of  cases. 
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Appendix  C:  The  TABS-MD  System 

TABS-MD  is  a  collection  of  generalized  computer  programs  and  utility 
codes  integrated  into  a  numerical  modeling  system.  TABS-MD  is  capable 
of  one-,  two-,  and/or  three-dimensional  computations,  but  only  the  one- 
and  two-dimensional  vertically  averaged  capability  will  be  discussed  in 
this  summary.  The  three-dimensional  version  of  the  code,  TABS-MDS 
(Multi-Dimensional  with  Sediment),  ERDC’s  version  of  RMAio,  will  not 
be  discussed  here.  The  system  is  used  for  studying  hydrodynamics,  sedi¬ 
mentation,  and  transport  problems  in  rivers,  reservoirs,  bays,  and  estu¬ 
aries.  A  schematic  representation  of  the  system  is  shown  in  Figure  Cl.  It 
can  be  used  either  as  a  stand-alone  solution  technique  or  as  a  step  in  a 
hybrid  modeling  approach.  The  basic  concept  is  to  calculate  water-surface 
elevations,  current  patterns,  sediment  erosion,  transport  and  deposition, 
the  resulting  bed  surface  elevations,  and  the  feedback  to  hydraulics.  Exist¬ 
ing  and  proposed  geometry  can  be  analyzed  to  determine  the  impact  on 
sedimentation  of  project  designs  and  to  determine  the  impact  of  project 
designs  on  salinity  and  on  the  stream  system.  The  system  is  described  in 
detail  by  Thomas  and  McAnally  (1985). 


The  three  basic  2D  depth-averaged  components  of  the  system  are  as 
follows: 

1.  “A  Two-Dimensional  Model  for  Free  Surface  Flows,”  RMA2. 

2.  “Sediment  Transport  in  Unsteady  2-Dimensional  Flows,  Horizontal 
Plane,”  SED2D. 

3.  “Two-Dimensional  Finite  Element  Program  for  Water  Quality,”  RMA4. 

The  U.S.  Army  Engineer  Research  and  Development  Center,  Coastal  and 
Hydraulics  Laboratory  (ERDC-CHL)  developed  and  maintains  TABS-MD. 
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RMA2  is  a  finite  element  solution  of  the  Reynolds  form  of  the  Navier- 
Stokes  equations  for  turbulent  flows.  Friction  is  calculated  with  Manning’s 
equation,  and  eddy  viscosity  coefficients  are  used  to  define  the  turbulent 
exchanges.  A  velocity  form  of  the  basic  equation  is  used  with  side  bounda¬ 
ries  treated  as  either  slip  or  static.  The  model  has  a  marsh  porosity  option 
as  well  as  the  ability  to  automatically  perform  wetting  and  drying.  Bound¬ 
ary  conditions  may  be  water-surface  elevations,  velocities,  discharges,  or 
tidal  radiation. 

The  sedimentation  model,  SED2D,  solves  the  convection-diffusion  equa¬ 
tion  with  bed  source-sink  terms.  These  terms  are  structured  for  either 
sand  or  cohesive  sediments.  The  Ackers  and  White  (1973)  procedure  is 
used  to  calculate  a  sediment  transport  potential  for  the  sands  from  which 
the  actual  transport  is  calculated  based  on  availability.  Clay  erosion  is 
based  on  work  by  Partheniades  (1962)  and  Ariathurai  et  al.  (1977)  and  the 
deposition  of  clay  used  Krone’s  equations.  Deposited  material  forms 
layers,  and  bookkeeping  allows  up  to  10  layers  at  each  node  for  maintain¬ 
ing  separate  material  types,  deposit  thickness,  and  age.  The  code  uses  the 
same  mesh  as  RMA2. 

Consistent  transport  calculations,  including  salinity,  are  made  under 
RMA4  using  a  form  of  the  convective-diffusion  equation  that  has  general 
source-sink  terms.  Up  to  six  conservative  substances  or  substances 
requiring  a  decay  term  can  be  routed.  The  code  uses  the  same  mesh  as 
RMA2.  The  model  accommodates  a  mixing  zone  outside  of  the  model 
boundaries  for  estimation  of  retrainment. 

Pre-  and  post-processing  and  analysis  of  TABS-MDS  models 

The  Surface  Water  Modeling  System  (SMS)  (BYU  2002)  is  a  comprehen¬ 
sive  graphical  user  environment  for  performing  model  conceptualizations, 
mesh  generation,  statistical  interpretation,  and  visual  examination  of  sur¬ 
face  water  model  simulations. 

SMS  is  a  pre-  and  post-processor  for  surface  water  modeling  and  analysis 
in  shallow  open  water  areas  such  as  rivers,  bays,  and  estuaries.  It  includes 
two-dimensional  finite  element,  two-dimensional  finite  difference,  three- 
dimensional  finite  element,  and  one-dimensional  step  backwater  modeling 
tools.  Interfaces  specifically  designed  to  facilitate  the  utilization  of  several 
numerical  models  comprise  the  modules  of  SMS.  Supported  models 
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include  the  TABS-MD  (GFGEN,  RMA2,  RMA4,  RMA10,  SED2D-WES) 
system. 

Each  TABS-MD  model  is  designed  to  address  a  specific  class  of  problem. 
RMA2  calculates  hydrodynamic  data  such  as  water  surface  elevations  and 
flow  velocities.  RMA4  tracks  contaminant  migration,  and  SED2D  calcu¬ 
lates  suspended  sediment  concentrations,  erosion,  and  deposition.  The 
models  support  both  steady-state  and  dynamic  analyses. 

The  finite  element  mesh  or  cross-section  entities,  along  with  associated 
boundary  conditions  necessary  for  analysis,  are  created  within  SMS  and 
then  saved  to  model-specific  files.  These  files  are  used  as  input  to  the 
hydrodynamic,  contaminant  migration,  and  sediment  transport  analysis 
engines.  The  numerical  models  create  solution  files  that  contain  the  water 
surface  elevations,  flow  velocities,  contaminant  concentrations,  sediment 
concentrations,  or  other  functional  data  at  each  node,  cell,  or  section. 

These  files  are  then  used  to  perform  the  analyses.  Resulting  solution  files 
can  be  read  into  SMS  to  generate  vector  plots,  color-shaded  contour  plots, 
time-history  diagrams,  and  solution  animation  sequences. 

Finite  element  modeling 

The  TABS-MD  numerical  models  employ  the  finite  element  method  to 
solve  the  governing  equations.  To  help  those  who  are  unfamiliar  with  the 
method  to  better  understand  the  system,  a  brief  description  of  the  method 
is  given  here. 

The  finite  element  method  approximates  a  solution  to  governing  equations 
by  dividing  the  area  of  interest  into  smaller  sub-areas,  which  are  called 
elements.  The  dependent  variables  (e.g.,  water-surface  elevations  or  sedi¬ 
ment  concentrations)  are  approximated  over  each  element  by  continuous 
functions  that  interpolate  based  on  unknown  point  (node)  values  of  the 
variables.  An  error,  defined  as  the  deviation  of  the  governing  equations 
using  the  approximate  solution  from  the  equation  using  the  correct  solu¬ 
tion,  is  minimized.  Then,  when  boundary  conditions  are  imposed,  a  set  of 
solvable  simultaneous  equations  is  created.  The  solution  is  continuous 
over  the  area  of  interest. 

In  one-dimensional  problems,  elements  are  line  segments.  In  two- 
dimensional  problems,  the  elements  are  polygons,  either  triangles  or 
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quadrilaterals.  Nodes  are  located  on  the  edges  of  elements  and  occasion¬ 
ally  inside  the  elements.  The  interpolating  functions  may  be  linear  or 
higher-order  polynomials.  Figure  C2  illustrates  a  quadrilateral  element 
with  eight  nodes  and  a  linear  solution  surface  where  F  is  the  interpolating 
function. 

Most  water  resource  applications  of  the  finite  element  method  use  the 
Galerkin  method  of  weighted  residuals  to  minimize  error.  In  this  method, 
the  residual— the  local  error  in  the  equations’  use  of  the  approximate 
solution— is  weighted  by  a  function  that  is  identical  to  the  interpolating 
function  and  then  minimized.  Minimization  results  in  a  set  of  simultane¬ 
ous  equations  in  terms  of  nodal  values  of  the  dependent  variable  (e.g., 
water-surface  elevations  or  sediment  concentration).  The  time  portion  of 
time-dependent  problems  can  be  solved  by  the  finite  element  method,  but 
it  is  generally  more  efficient  to  express  derivatives  with  respect  to  time  in 
finite  difference  form. 
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Hydrodynamic  model,  RMA2 

Applications 

RMA2  is  designed  for  far-field  problems  in  which  vertical  accelerations  are 
negligible  and  the  velocity  vectors  at  a  node  generally  point  in  the  same 
directions  over  the  entire  depth  of  the  water  column  at  any  instant  of  time. 
It  expects  a  vertically  homogeneous  fluid  with  a  free  surface.  The  model 
will  define  the  response  to  a  specified  horizontally  inhomogeneous  fluid. 
Both  steady  and  unsteady  state  problems  can  be  analyzed.  A  surface  wind 
stress  can  be  imposed  and  the  effects  of  the  earth’s  rotation  (Coriolis 
effect)  can  be  included. 


RMA2  has  been  applied  to  calculate  water  levels  and  flow  distribution 
around  islands;  flow  at  bridges  having  one  or  more  relief  openings,  in  con¬ 
tracting  and  expanding  reaches,  into  and  out  of  off-channel  hydropower 
plants,  at  river  junctions,  and  into  and  out  of  pumping  plant  channels; 
circulation  and  transport  in  water-bodies  with  wetlands;  and  general  water 
levels  and  flow  patterns  in  rivers,  reservoirs,  and  estuaries. 

Limitations 

RMA2  is  not  designed  for  near-field  problems  where  flow  structure  inter¬ 
actions  (such  as  vortices,  vibrations,  or  vertical  accelerations)  are  of 
interest.  Areas  of  vertically  stratified  flow  are  beyond  this  program’s 
capability  unless  it  is  used  in  a  hybrid  modeling  approach.  It  is  two- 
dimensional  in  the  horizontal  plane,  so  zones  where  the  bottom  current  is 
in  a  different  direction  from  the  surface  current  must  be  analyzed  with 
care.  It  is  a  free-surface  calculation  for  sub-critical  flow  problems. 


Governing  equations 

The  generalized  computer  program  RMA2  solves  the  depth-integrated 
equations  of  fluid  mass  and  momentum  conservation  in  two  horizontal 
directions.  The  form  of  the  solved  equations  is 


,  du  ,  du  ,  du  h 

h - b  hu - b  hv - 

dt  dx  dy  p 


d2u  d2u 


V 


+  gh 


da  dh 
dx  dx 


H - - 5-  (u2  -b  v2)  1  -  Z,  V2  cos  \p  -  2/200  u  sin  cp  =  0 

(1.486/r1/6)2  1  j 


(Cl) 
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h  =  depth 

ii, v  =  x  and  y  direction  velocities,  respectively 
x,y,t  =  Cartesian  coordinates  and  time 
p  =  density  of  fluid 
8  =  eddy  viscosity  coefficient, 

for  xx  =  normal  direction  on  x-axis  surface; 
yy  =  normal  direction  on  y-axis  surface; 
xy  and  yx  =  shear  direction  on  each  surface 
g  =  acceleration  due  to  gravity 
a  =  elevation  of  bottom 
n  =  Manning’s  n  value 

1.486  =  conversion  from  SI  (metric)  to  non-SI  units 
£  =  empirical  wind  shear  coefficient 
Va  =  wind  speed 
\p  =  wind  direction 
co  =  rate  of  earth’s  angular  rotation 


cp  =  local  latitude. 


Equations  Cl,  C2,  and  C3  are  solved  by  the  finite  element  method  using 
Galerkin  weighted  residuals.  The  elements  may  be  one-dimensional  lines 
or  two-dimensional  quadrilaterals  or  triangles  and  may  have  curved  (para¬ 
bolic)  sides.  The  shape  functions  are  quadratic  for  velocity  and  linear  for 
depth.  Integration  in  space  is  performed  by  Gaussian  integration.  Deriva¬ 
tives  in  time  are  replaced  by  a  nonlinear  finite  difference  approximation. 
Variables  are  assumed  to  vary  over  each  time  interval  in  the  form 


f(t)=  /( 0)  +  at  +  btc  t0<  t  Z  t0  +  A t 


(C4) 
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which  is  differentiated  with  respect  to  time,  and  cast  in  finite  difference 
form.  Letters  a,  b,  and  c  are  constants.  It  has  been  found  by  experiment 
that  the  best  value  for  c  is  1.5  (Norton  and  King  19 77). 

The  solution  is  fully  implicit  and  the  set  of  simultaneous  equations  is 
solved  by  Newton-Raphson  nonlinear  iteration. 

Sediment  transport  model,  SED2D 

Applications 

SED2D  can  be  applied  to  clay  and/or  sand  bed  sediments  where  flow 
velocities  can  be  considered  two-dimensional  (i.e.,  the  speed  and  direction 
can  be  satisfactorily  represented  as  a  depth-averaged  velocity).  It  is  useful 
for  both  deposition  and  erosion  studies  and,  to  a  limited  extent,  for  stream 
width  studies.  The  program  treats  two  categories  of  sediment:  noncohe- 
sive,  which  is  referred  to  here  as  sand,  and  cohesive,  which  is  referred  to  as 
clay. 

Limitations 

Both  clay  and  sand  may  be  analyzed,  but  SED2D  considers  a  single,  effec¬ 
tive  grain  size  for  each  and  treats  each  separately.  Fall  velocity  must  be 
prescribed  along  with  the  water-surface  elevations,  x-velocity,  y-velocity, 
diffusion  coefficients  and,  for  clay  sediment,  bed  density,  critical  shear 
stresses  for  erosion,  erosion  rate  constants,  and  critical  shear  stress  for 
deposition. 

The  program  does  not  compute  water-surface  elevations  or  velocities,  so 
those  data  must  be  provided.  For  complicated  geometries,  the  numerical 
model  for  hydrodynamic  computations  (i.e.,  RMA2)  is  used.  However,  at 
this  time,  SED2D  can  accept  only  a  two-dimensional  network. 

Governing  equations 

SED2D  solves  the  depth-integrated  convection-dispersion  equation  in  two 
horizontal  dimensions  for  a  single  sediment  constituent.  For  a  more  com¬ 
plete  description,  see  Appendix  G  of  Thomas  and  McAnally  (1985).  The 
form  of  the  solved  equation  is 
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dc  dc  dc 

- h  u - b  v — 

dt  dx  dy 


d_ 

dx 


D. 


dC] 

dx 


d_ 
+  dy 


D, 


dc_ 

dy 


-\-  cqC  T-  0-2  —  0 


(C5) 


where: 

C  =  concentration  of  sediment 
u  =  depth-integrated  velocity  in  x-direction 
Dx  =  dispersion  coefficient  in  x-direction 
v  =  Depth  integrated  velocity  in  y-direction 
Dy  =  dispersion  coefficient  in  y-direction 
di  =  coefficient  of  concentration-dependent  source/sink  term 
a2  =  coefficient  of  source/sink  term. 


The  source/sink  terms  in  Equation  C5  are  computed  in  routines  that  treat 
the  interaction  of  the  flow  and  the  bed.  Separate  sections  of  the  code 
handle  computations  for  clay  bed  and  sand  bed  problems. 

Sand  transport 

The  source/sink  terms  are  evaluated  by  first  computing  a  potential  sand 
transport  capacity  for  the  specified  flow  conditions,  comparing  that 
capacity  with  the  amount  of  sand  actually  being  transported,  and  then 
eroding  from  or  depositing  to  the  bed  at  a  rate  that  would  approach  the 
equilibrium  value  after  sufficient  elapsed  time. 


The  potential  sand  transport  capacity  in  the  model  is  computed  by  the 
method  of  Ackers  and  White  (1973),  which  uses  a  transport  power  (work 
rate)  approach.  It  has  been  shown  to  provide  superior  results  for  transport 
under  steady-flow  conditions  (White  et  al.  1975)  and  for  combined  waves 
and  currents  (Swart  1976).  Flume  tests  at  ERDC  have  shown  that  the 
concept  is  valid  for  transport  by  estuarine  currents. 

The  total  load  transport  function  of  Ackers  and  White  is  based  upon  a 
dimensionless  grain  size 


D 


g(s-i) 
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1/3 
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where: 

D  =  sediment  particle  diameter 
s  =  specific  gravity  of  the  sediment 
v  =  kinematic  viscosity  of  the  fluid 

and  a  sediment  mobility  parameter 


F*r  = 


(1  ~n) 


pgD(s-l) 


1/2 


(C7) 


where: 

t'  =  total  boundary  shear  stress  =  p gRS 
R  =  hydraulic  radius 
S  =  slope  of  water  surface 

n  =  a  coefficient  expressing  the  relative  importance  of  bed-load 
and  suspended-load  transport,  given  in  Equation  C9. 

NOTE: 

n  =  1  for  fine  sediments 
n  =  o  for  coarse  sediments 
t  =  boundary  surface  shear  stress. 

The  surface  shear  stress  is  that  part  of  total  shear  stress  attributable  to  the 
rough  surface  of  the  bed  only,  i.e.,  not  including  the  part  attributable  to 
bed  forms  and  geometry.  It  therefore  corresponds  to  the  shear  stress  that 
the  flow  would  exert  on  a  plane  bed. 

The  total  sediment  transport  is  (in  kg/m3)  expressed  as  an  effective 
concentration 


m 

^ G 

Pf7 

h 

X 

(C8) 


where  U  is  the  average  flow  speed,  and  for  1  <  Dgr  <60 
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log  Ca 


n  =  1.00  -  0.56  log Dgr 

(C9) 

A  -  °^1+0.14 

Hr 

(CIO) 

2.86  log  Dgr-  (log  Dgrf~  3.53 

(Cll) 

9  66 

m  =  +  1.34 

Dgr 

(C12) 

For  Dgr  <  6o 


n  =  0.00  (C13) 

A  =  0.17  (C14) 

Ca  =  0.025  (C15) 

m  =  1.5  (C16) 


Note  the  Ca  has  units  consistent  with  Gp  (kg/m3  for  SED2D). 


Equations  C6  through  C16  result  in  a  potential  sediment  concentration  Gp. 
This  value  is  the  depth-averaged  concentration  of  sediment  that  will  occur 
if  an  equilibrium  transport  rate  is  reached  with  a  limited  supply  of  sedi¬ 
ment.  The  rate  of  sediment  deposition  (or  erosion)  is  then  computed  as 


R 


(C 17) 


where: 

C  =  present  sediment  concentration 
tc  =  time  constant. 


ERDC/CHL  TR-08-11 


94 


For  deposition,  the  time  constant  is 


t 


C 


At 


larger  of 


or 


C^h 

Vs 


(C18) 


and  for  erosion  it  is 
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C 


At 


larger  of 


or 


Cjl_ 

U 


(C19) 


where: 

At  =  computational  timestep 

Cd  =  response  time  coefficient  for  deposition 

Vs  =  sediment  settling  velocity 

Ce  =  response  time  coefficient  for  erosion. 

The  sand  bed  has  a  specified  initial  thickness  which  limits  the  amount  of 
erosion  to  that  thickness. 

Cohesive  sediments  transport 

Cohesive  sediments  (usually  clays  and  some  silts)  are  considered  to  be 
depositional  if  the  bed  shear  stress  exerted  by  the  flow  is  less  than  a  critical 
value  t d.  When  that  value  occurs,  the  deposition  rate  is  given  by  Krone’s 
(1962)  equation 


S 


forC  <  Cc 
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where: 

S  =  source  term 

Vs  =  fall  velocity  of  a  sediment  particle 
h  =  flow  depth 

C  =  sediment  concentration  in  water  column 
x  =  bed  shear  stress 
t d  =  critical  shear  stress  for  deposition 
Cc  =  critical  concentration  =  300  mg/L. 

If  the  bed  shear  stress  is  greater  than  the  critical  value  for  particle  erosion 
Te,  material  is  removed  from  the  bed.  The  source  term  is  then  computed  by 
Ariathurai’s  (Ariathurai  et  al.  1977)  adaptation  of  Partheniades’  (1962) 
findings: 


S 


P 

h 


^  -1 


for  x  >  xe 


(C21) 


where  P  is  the  erosion  rate  constant,  unless  the  shear  stress  is  also  greater 
than  the  critical  value  for  mass  erosion.  When  this  value  is  exceeded,  mass 
failure  of  a  sediment  layer  occurs  and 


S  =  for  x  >  t 
hAt 


(C22) 


where: 

Tl  =  thickness  of  the  failed  layer 
p l  =  density  of  the  failed  layer 
At  =  time  interval  over  which  failure  occurs 
t s  =  bulk  shear  strength  of  the  layer. 

The  cohesive  sediment  bed  consists  of  1  to  10  layers,  each  with  a  distinct 
density  and  erosion  resistance.  The  layers  consolidate  with  overburden 
and  time. 

Bed  shear  stress 

Bed  shear  stresses  are  calculated  from  the  flow  speed  according  to  one  of 
four  optional  equations:  the  smooth-wall  log  velocity  profile  or  Manning 
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equation  for  flows  alone;  and  a  smooth  bed  or  rippled  bed  equation  for 
combined  currents  and  wind  waves.  Shear  stresses  are  calculated  using  the 
shear  velocity  concept 


Tb  =  pul 


(C23) 


where: 

Tb  =  bed  shear  stress 
u»  =  shear  velocity 


and  the  shear  velocity  is  calculated  by  one  of  four  methods: 


l.  Smooth-wall  log  velocity  profiles 


u 


u» 


5.75  log 


3.23 


u*h 


(C24) 


Equation  C25  is  applicable  to  the  lower  15  percent  of  the  boundary 
layer  when 


—  >30  (C25) 

v 


where  u  is  the  mean  flow  velocity  (resultant  of  u  and  v  velocity 
components). 


2.  The  Manning  shear  stress  equation 

_  (hn)V^ 

t-t*  -  't  /r 

CME  (h)1 


(C26) 


where  CME  is  a  coefficient  of  1  for  SI  (metric)  units  and  1.486  for 
English  units  of  measurement,  and  n  is  the  Mannings  n-value. 

3.  A  Jonsson-type  equation  for  surface  shear  stress  (plane  beds)  caused  by 
waves  and  currents 
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U»  = 


fwUom  +  fc11 


Uom  +U 


u+uomy 


(C27) 


where: 

fw  =  shear  stress  coefficient  for  waves 
Uom  =  maximum  orbital  velocity  of  waves 
fc  =  shear  stress  coefficient  for  currents. 

4.  A  Bijker-type  equation  for  total  shear  stress  caused  by  waves  and  current 


U, 


r  — 2  1  r 

fcu  +4  fwU 
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(C28) 


Solution  method 

Equation  C5  is  solved  by  the  finite  element  method  using  Galerkin 
weighted  residuals.  Like  RMA2,  which  uses  the  same  general  solution 
technique,  elements  are  quadrilateral  or  triangular  and  may  have  para¬ 
bolic  sides.  Shape  functions  are  quadratic.  Integration  in  space  is 
Gaussian.  Time-stepping  is  performed  by  a  Crank-Nicholson  approach 
with  a  weighting  factor  (0)  of  0.66.  A  front-type  solver  similar  to  that  in 
RMA2  is  used  to  solve  the  simultaneous  equations. 

Water  quality  transport  model,  RMA4 

Applications 

The  water  quality  transport  model,  RMA4,  is  designed  to  simulate  the 
depth-average  advection-diffusion  process  in  most  water  bodies  with  a 
free  surface.  The  model  is  used  for  investigating  the  physical  processes  of 
migration  and  mixing  of  a  soluble  substance  in  reservoirs,  rivers,  bays, 
estuaries,  and  coastal  zones.  The  model  is  useful  for  evaluation  of  the  basic 
processes  or  for  defining  the  effectiveness  of  remedial  measures.  For 
complex  geometries  the  model  utilizes  the  depth-averaged  hydrodynamics 
from  RMA2. 

The  water  quality  model  has  been  applied  to  define  the  horizontal  salinity 
distribution,  trace  temperature  effects  from  power  plants,  calculate  resi¬ 
dence  times  of  harbors  or  basins,  optimize  the  placement  of  outfalls, 
identify  potential  critical  areas  for  oil  spills  or  other  pollutants  spread, 
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evaluate  turbidity  plume  extent,  and  monitor  other  water  quality  criterion 
within  game  and  fish  habitats. 

Limitations 

The  formulation  of  RMA4  is  limited  to  one-dimensional  (cross-sectionally 
averaged)  and  two-dimensional  (depth-averaged)  situations  in  which  the 
concentration  is  fairly  well  mixed  in  the  vertical.  It  will  not  provide  accu¬ 
rate  concentrations  for  stratified  situations  in  which  the  constituent  con¬ 
centration  influences  the  density  of  the  fluid.  In  addition,  the  accuracy  of 
the  transport  model  is  dependent  on  the  accuracy  of  the  hydrodynamics 
(e.g.,  as  supplied  from  RMA2  or  another  flow  solution). 


Governing  equations 

The  ERDC  version  of  RMA4  is  a  revised  version  of  RMA4  as  developed  by 
King  and  Rachide  (1989).  The  generalized  computer  program  solves  the 
depth-integrated  equations  of  the  transport  and  mixing  process.  The  form 
of  the  equations  solved  is: 


dc  dc  dc  d  dc  d  dc  , 

dt  dx  dy  dx  dx  dy  y  dy 


=  0  (C29) 


where: 
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=  water  depth 
=  constituent  concentration 
=  time 

=  velocity  components 
=  turbulent  mixing  coefficients 
=  first  order  decay 
=  source/sink  of  constituent. 


Note  that  the  basic  governing  equation  for  RMA4  is  the  same  as  for  the 
sediment  transport  model,  SED2D.  The  difference  between  the  two  models 
lies  in  the  source/sink  terms. 


Equation  C29  is  solved  by  the  finite  element  method  using  Galerkin 
weighted  residuals.  As  with  the  hydrodynamic  model,  RMA2,  the  trans¬ 
port  model  RMA4  handles  one-dimensional  segments  or  two-dimensional 
quadrilaterals  or  triangles  with  the  option  for  curved  sides.  Spatial 
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integration  of  the  equation  is  performed  by  Gaussian  techniques  and  the 
temporal  variations  are  handled  by  nonlinear  finite  differences,  consistent 
with  the  method  described  for  RMA2. 

The  boundary  conditions  for  RMA4  are  specified  in  several  optional  ways. 
The  boundary  concentration  may  be  specified  absolutely  at  a  certain  level 
regardless  of  the  flow  direction;  the  concentration  can  be  specified  to  be 
applied  only  when  the  water  is  leaving  the  model;  or  a  mixing  zone  may  be 
specified  just  beyond  the  model  boundary  to  provide  the  possibility  of  re¬ 
entrainment  of  constituent  into  the  model  that  may  have  crossed  the 
boundary  earlier.  For  a  more  detailed  description  of  the  constituent  trans¬ 
port  model,  RMA4,  see  King  and  Rachiele  (1989). 

Within  the  one-dimensional  formulation  of  the  model,  there  is  a  provision 
for  defining  the  constituent  concentration  mixing  and  transport  at  control 
structures  as  they  may  have  been  specified  in  RMA2.  This  allows  for  either 
a  flow-through  condition,  as  for  example  for  a  weir-type  flow,  or  for  a 
mixing  chamber  type  of  flux,  which  would  be  appropriate  for  a  navigation 
lock. 
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